Calculating spectral signatures for each region

AFNI Message Board
1

Hello,

I have been trying to calculate spectral features of resting-state scans for low frequency (0.01-0.1Hz) and high frequency (0.1-0.5 Hz) bands using Matlab but wanted to see if there might be a set of AFNI commands that may be more convenient.

My scans consist of 385 volumes with a TR of 1 second so Nyquist frequency would be 0.5 Hz. This older post on AFNI discussion board was super helpful and wanted to see if the 'pipeline' below looks okay (3dPeriodogram and 1dFFT - How exactly is the power computed?).

  • Use 3dPeriodogram on every voxel in AFNI on the preprocessed time-series, using the command below:
    3dPeriodogram
    -prefix p
    -taper 0.1
    -nfft 512
    rest.nii.gz

  • Cut the 3dPeriodogram results to the desired frequency range via 3dTcat, using commands below.
    3dTcat -prefix high_freq p+orig'[51..255]'
    3dTcat -prefix low_freq p+orig'[0..50]'

  • Do voxel-wise mean respectively for the low- and high frequency files.

3dmaskave -mask mask low_freq > roi_mean_lowfreq.1D

  • And then if I want to get an average power value for 0.01-0.1Hz band then I would just average all values for a given ROI from the 1D file I get from 3dmaskave.

Just wanted to see if I could get any input on whether this makes sense.

Thank you in advance.

2

Howdy-

If you are interested in RSFC parameters, there are programs specifically for this that might be useful. Getting the average power value is kind of like the ALFF (amplitude of low frequency fluctuations, which is a sum of amplitudes across the specified band range) parameter (power is the square of amplitude), and you could get that for specific bands. The fALFF parameter (fractional ALFF) would give you the ratio of ALFF in your specified range to the sum of all amplitudes in the time series.

One big question will be: did you include censoring in your processing?

  • If you did not, you can use 3dRSFC to estimate a set of RSFC parameters.
  • But if censoring was included, then there are "gaps" in your data, and these must be navigated when using an spectrum-based parameter, like ALFF and fALFF. When censoring has been used, the standard Fourier Transform does not apply, because it assumes uniform sampling in the input time series; censoring destroys that uniformity. Instead, one needs to use the Lomb-Scargle Transform to estimate a frequency spectrum and then get RSFC parameters from those results. In AFNI, this is done with a combination of 3dLombScargle and 3dAmpToRSFC. These programs are introduced+demonstrated a bit in this ISMRM-2018 poster:
    Lomb-Scargle your way to RSFC parameter estimation in AFNI-FATCAT

You could use the above programs on volumetric datasets to get voxelwise estimates. Note that another consideration is that as more and more censoring occurs, there will be a bit of a bias potential even within the Lomb-Scargle Transform---see the above poster. If censoring is more than, like 10% of your time series points, that might be a notable issue.

If you need this ROI-based measures, I guess you would first calculate the average time series of each ROI, and store those in a 1D file, and then I think you can put that into either 3dRSFC or 3dLombScargle+3dAmpToRSFC. In practice, getting voxelwise estimates of your RSFC parameters per voxel and then averaging those within an ROI might yield very similar results. I am not sure if there is a deep theoretical reason to average the time series and then calculate a parameter, or calculate the parameters and average those results.

--pt

3

Hi Dr. Taylor,

Thank you for the suggestions! I do not censor individual bricks (but instead simply exclude participants with too much overall movement), so 3dRSFC should be applicable.
But, I just wanted to go over a few things so that I understand the issue more clearly. Hope this is okay.

  1. What's the fundamental difference between what I proposed to do using 3dPeriodogram vs. 3dRSFC as you suggested? My understanding is the former gives you power by fine-grained frequency while the latter you can get amplitude values for pre-specified frequency range.

  2. 3dRSFC seems to rely on 3dBandPass, which I prefer not to use because it would automatically detrend the data (the data's already de-trended). Would there be a way to skip detrending and also use 3dRSFC?

  3. I want to draw power spectrograms and spectral density plots by region like the one here in Figure 1C. I don't thin 3dRSFC would allow me to do that... if the steps I described in my initial post seem incorrect, would you be able to let me know what you think would be the best approach?

  4. I am also hoping to create a 3d map of power values so that I can do voxel-wise analyses, so that I can create average maps by group. Would the steps I had outlined in my initial post allow me to do this?

Thank you again. Happy new year.

4

Howdy-

OK, that's all good to know. A key branching point for decision making here is whether or not censoring has been done (I'm just restating this for future generations and/or chatbot scrapers):

  • if censoring has been done, use 3dLombScargle to get the frequency spectrum estimate (via the eponymous transform), and then 3dAmpToRSFC for the RSFC params from that.
  • if censoring has not been done, we are in Fourier Transform-able territory, and 3dPeriodogram and 3dRSFC are both good to go.

So, then to your points/questions:

  1. Sure, using 3dPeriodogram is fine. It basically does an intermediate step done by 3dRSFC, giving you the frequency spectrum. If that be of interest, go ahead and use 3dPeriodogram, sure.

  2. 3dRSFC has a -nodetrend option to skip detrending. As a historical note, 3dRSFC was built by starting with a copy of 3dBandpass and then adjusting it, adding/subtracting options, etc. So, it doesn't call it under the hood, it just has very similar style/usage.

  3. Use can use 3dPeriodogram, sure. The output of that has a "frequency series" at each voxel, and the delta_frequency is given by "3dinfo -tr .."; if you look at the 3dinfo output on that dataset, you will see that the units of the "Time step" are in fact Hz, as appropriate for frequencies.

  4. 3dPeriodogram will give you a 4D dataset, with a frequency series at each voxel. There are many ways to collapse/calculate values across the frequency axis. Essentially, RSFC parameters are one particular set. But you can use 3dTstat for many others, too.

How does that seem?

And happy end of 2025/start of 2026 to you, as well!

--pt

5

Dear Dr. Taylor,

I tried a few things with 3dRSFC and wanted to see if I could get your advice on understanding different pieces.

  1. I tried below both without and with the -mask option, and it seemed to yield the same voxelwise values. Also, even when I included the -mask option, the log still reported ++ No mask ==> processing all 2097152 voxels . In addition, although I included the -nodetrend option, the log still says + -- quadratic detrend . Could I get your thoughts on what exactly the -mask and -nodetrend options are doing here? Is it possible that I am using this command in the wrong way? I have attached the log below.
3dRSFC -prefix rsfc 0.01 0.1 pb03.id.r01.blur.nii -nodetrend -mask follow_ROI_gm.nii

++ 3dRSFC (from 3dBandpass by RW Cox): version THETA: AFNI version=AFNI_24.0.12 (Mar 12 2024) [64-bit]
++ Authored by: PA Taylor
++ Data length = 385 FFT length = 386
+ bandpass: ntime=385 nFFT=386 dt=1 dFreq=0.00259067 Nyquist=0.5 passband indexes=4..39
++ Loading input dataset time series
++ No mask ==> processing all 2097152 voxels
++ Checking dataset for initial transients [use '-notrans' to skip this test]
+ No widespread initial positive transient detected :-)
++ THD_dset_to_vectim: allocated 3229614080 bytes
++ Bandpassing data time series
+ -- quadratic detrend
+ -- FFTs
++ Creating output dataset in memory, then writing it
++ Output dataset ./rsfc_LFF+orig.BRIK
++ THD_dset_to_vectim: allocated 3229614080 bytes
++ Bandpassing data time series
+ -- quadratic detrend
++ Creating output dataset 2 in memory
++ Starting the (f)ALaFFel calcs
  1. I also tried analyzing higher frequencies by writing: 3dRSFC -prefix highfreq 0.11 0.44 pb03.id.r01.blur.nii -nodetrend -mask follow_ROI_gm.nii . If I want to examine amplitude specific to this frequency range, would it be correct to focus on the RSFA output files? In addition, this command still generated ALFF and LFF files, which is confusing to me because my understanding is that these measures typically correspond to the 0.01 to 0.08 Hz band, which lies outside the frequency range I specified. I would appreciate any clarification on how these outputs are being computed in this context.

Thank you so much!

6

Hi-

Re. Q1, about the mask:
This program has usage style:

3dRSFC [options] fbot ftop dataset

So, the input dataset has to be the last thing on the command line, and it must be preceded by 2 numbers that define the bandpassing range. All options (like -mask .. and nodetrend) have to be placed to the left of/before those 3 args.

Your command was written as:

3dRSFC -prefix rsfc 0.01 0.1 pb03.id.r01.blur.nii -nodetrend -mask follow_ROI_gm.nii

... which doesn't do that. So, it was good that you checked the output text. Please try writing it like:

3dRSFC -prefix rsfc -nodetrend -mask follow_ROI_gm.nii 0.01 0.1 pb03.id.r01.blur.nii

Re. Q2, about higher frequencies:
You can certainly check higher frequency ranges, up to the Nyquist frequency, which is defined in terms of the sampling time (= TR here) as:
f_nyquist = 1/(2*TR)
So, for a TR = 2s FMRI acquisition, f_nyquist = 0.25 Hz.

You might be interested in this nice paper, that looked at information across different frequency bands, including both within the standard LFF range (~0.001 - 0.1 Hz) and above:

  • Gohel SR, Biswal BB. Functional integration between brain regions at rest occurs in multiple-frequency bands. Brain Connect. 2015 Feb;5(1):23-34. doi: 10.1089/brain.2013.0210.

Whatever frequency range you pick, you will get things called ALFF and LFF. It's true those might not correspond to exactly what is generally considered a "low frequency range", they just borrow/follow the same mathematical formulations how how ALFF and LFF are calculated, for the range the user provides. (If it's helpful, you can see a set of the formulations in Appendix A of this paper, which includes RSFA, which is very similar to ALFF; RSFA is basically just an L2 norm while ALFF is an L1 norm, of the same input.) So, if you use a frequency range of 0.11-0.44 Hz, the "ALFF" parameter will be the (scaled) sum of amplitudes of frequency magnitudes within that range. Basically, the quantities/formulations are the same as what the standard usage of terms is, but if users redefine what the frequency range is, that's fine.

How does that all seem?

--pt

7

Dear Dr. Taylor,

These all seem great. 3dRSFC worked for me. And thank you for clarifying the difference between ALFF and RSFA.

  1. I also tried calculating power for low (0.01-0.1 Hz) and higher (0.11-0.44 Hz) frequency bands using the code below. Would you be able to let me know if they look okay to you? My knowledge is that what 3dPeriodogram outputs in each brick is power spectral density (squared amplitude of the Fourier transform at each respective frequency, normalized by the FFT length), and I've written codes below based on this understanding. errs.id.tproject+orig is I believe the de-trended output of proc_afni.py.
3dPeriodogram -prefix pdgram -nfft 512 errts.id.tproject+orig
3dTstat -mean -prefix avg_lowfreq.nii.gz 'pdgram+orig[5..51]'
3dTstat -mean -prefix avg_highfreq.nii.gz 'pdgram+orig[56..224]'

##below is just outputting files
echo -n "${subject} " >> $outgm
3dmaskave -mask follow_ROI_gm.nii avg_lowfreq.nii.gz >> $outgm
echo -n "${subject} " >> $outroi
3dROIstats -mask follow_ROI_aparcaseg.nii avg_lowfreq.nii.gz >> $outroi

echo -n "${subject} " >> $outgmh
3dmaskave -mask follow_ROI_gm.nii avg_highfreq.nii.gz >> $outgmh
echo -n "${subject}" >> $outroih
3dROIstats -mask follow_ROI_aparcaseg.nii avg_highfreq.nii.gz >> $outroih
  1. I also want to do some voxelwise group comparisons of the low and high frequency power brain images. After noticing that AFNI files take up less space, I've been trying to use AFNI instead of fsl-flirt to warp EPI to MNI space. Would the codes I used below seem ok?

The part I'm a bit confused about is EPI resolution is 2x2x2mm while ANAT is 1x1x1. I want to use warp generated by aligning ANAT to MNI to ultimately align EPI to MNI, but still want to keep the MNI-aligned EPI in 2x2x2mm resolution.

@SSwarper                                                       \
  -input  anat                                        \
  -base   MNI152_2009_template_SSW.nii.gz  \                                                  \
  -subid  participant

cat_matvec -ONELINE orig_al_junk_mat.aff12.1D -I > epi2anat.aff12.1D

3dNwarpCat \
  -prefix epi2mni_WARP.nii.gz \
  -warp1  SSW_WARP.nii.gz \
  -warp2  SSW_aff12.1D \
  -warp3  epi2anat.aff12.1D

3dresample -input $mni -dxyz 2 2 2 -prefix mni_2mm.nii.gz

3dNwarpApply \
  -source $epi_3d \
  -master mni_2mm.nii.gz \
  -nwarp  epi2mni_WARP.nii.gz \
  -interp wsinc5 \
  -prefix epi_MNI_2mm.nii.gz
8

Howdy-

There are a few topics here.

Re. 1/3dPeriodogram: What is your TR? The upper frequency possible to reconstruct is determined by the Nyquist frequency, which I mentioned above. On the periodogram side, the delta_f (frequency step size between successive volumes in the periodogram) will be given by the total time, start-to-finish, of the input time series:
delta_f = TR * N_timepoints
Indeed, you can zeropad the time series to increase N_timepoints; this adjusts the sampling rate of the frequency spectrum by shrinking delta_f. Given your stated frequency bands and subbrick selectors, I guess you have TR~1s, so you should be OK for those ranges?

Also, I guess there is a gap in the subbrick selectors between the low and higher bands (between index 51 and index 56) to leave some room for taper dropoff?

Your calculations with 3dTstat calculate is the total power in a band divided by the number of time points in the band---so, average power per band. That is fine to also average within each ROIs or GM mask, sure.

Re. 2/warping: I'm a bit confused as to what the datasets are. By seeing the name errts* as the input to Q1's 3dPeriodogram command, I assumed this was your final processed data from afni_proc.py, and that as such this would be in the final desired space. But the mention of nonlinear warping after having the errts* (which has a "+orig" in its file name here, so I can see it is in original/subject space) makes it seem like that is not the case.

I would like to suggestion running nonlinear warping before running afni_proc.py (you can use @SSwarper for this, or what I would recommend instead is the updated version sswarper2, which has the same set of outputs and usage), and then providing those results to afni_proc.py, to include the nonlinear warping to standard space within that full subject processing pipeline. This reduces the total number of regriddings that the EPI data goes through, so that extra+unnecessary blurring is not incurred beyond what is specified in processing. That should also help give you all your data on the same grid, across the group, esp. in you include -volreg_warp_dxyz VALUE to explicitly choose a final EPI voxel size in the output space. The work you are doing separately with 3dNwarpCat and 3dNwarpApply would all be handled more simply, directly and streamlinedly within afni_proc.py. Several afni_proc.py help examples demo doing this, as it is our recommended way to go for volumetric analyses in standard space.

Note that in terms of choosing output voxel size when using 3dNwarpApply, there is the -dxyz .. option, so you don't have to regrid the master and use that to specify voxel dims (although that can be done).

In terms of using 3dNwarpCat with multiple warps and sswarper output, that is the correct set of warps from sswarper, and also if you have been using afni_proc.py for the initial processing and using options I would guess are recommended, that guessing from those file names, that seems appropriate. But again, it would be better to let afni_proc.py manage all of this internally and in a single step.

--pt