what measure is recommended to summarize "spikiness" in an area

[carolin31]

Hello afni experts,

I am focusing in my analysis on the HIP and I would like to assess the “spikiness” in the left and right Hippocampus as a quality check in order to exclude subjects that have too many spikes in the HIP. I used 3dDespike with the -ssave option to get the spikiness measure s for each voxel. This is a 3D+time data set. What measure would you recommend for me to extract to summarize the “spikiness” in all voxels of the HIP in one measure over time?
Is there a way to count the number of spikes? Or to calculate some kind of mean deviation measure from the mean signal?

Thank you very much for any help
Carolin

[ptaylor]

Hi, Carolin-

With standard FMRI processing with afni_proc.py, one criterion to censor a particular time point is to:

• calculate outliers at each voxel, across all time, within the brain mask
• then count the fraction of outliers within the brain mask (and compare that fractional value to a user-chosen threshold, to decide about censoring there).

That sounds like what you want to do, per ROI? You could use either outlier counts per ROI, or outlier fraction per ROI, per time point?

The following tcsh script loops over two ROIs in an “ROI dset”, and calculates the outlier fraction within the “EPI dset” at each point in time, and it dumps that result to a 1D/time-series file; it also simply plots each:

#!/bin/tcsh

set dset_roi = # user puts a dset name here
set dset_epi = # user puts a dset name here

# which ROIs to select (by ROI value)
set all_roival = ( 1 2 )

# loop over each ROI value
foreach roival ( ${all_roival} )
    # compute outlier fraction for each volume, 
    # selecting an ROI value from the $dset_roi to be a mask 
    3dToutcount                                           \
        -mask ${dset_roi}"<${roival}>"                    \
        -fraction -polort 1 -legendre                     \
        ${dset_epi} > out_frac_${roival}.txt
    # plot
    1dplot out_frac_${roival}.txt &
end

You can remove “-fraction” within 3dToutcount, if you want a number of outliers per time point. I “stole” the command from an afni_proc.py-generated script, using the polort level and legendre polynomial aspects.

–pt

[carolin31]

Thank you!
I was able to successfully create a 1D file for each subject that shows the outlier fraction for each time point for the left and then the right HIP.

How can I now calculate the number of timepoints/volumes that exceed a certain threshold for outliers in the 1D file?
For my analysis I need to have a temporal continuity of the data. This means I can not censor individual volumes. Instead I want to use the number of volumes that would need to be censored based on outlier fraction as a criteria for whether to exclude a subject or not.

Also how is the count calculated? Is it the number of spikes averaged across all voxels in my ROI for each timepoints? If so isn t that similar to the fraction? In what instances is it recommended to use the counts and in what the fraction?

Carolin

[carolin31]

Hello AFNI experts,

I was writing before to ask what metric I could use to summarize the spikiness measure in an ROI after using 3dDespike.

What I did so far is using 3dDespike which gave me an output dataset that includes the spikiness measure ‘s’ for each timepoints of every voxel.
I then used 3dToutcount to calculate the fraction of outliers at each time point in my ROI.

However, I realized just now that this might not be the right function to use as 3dToutcount calculates a trend and MAD across the time series for each voxel and then determines the outliers when a time point is far away from this trend. So in my case it calculates in each voxel a trend over the spikiness measure ‘s’ across all time points and outlier voxels will be timepoints where the spikiness measure deviates from the trend in this voxel.That does not sound like the right thing to do.

What I would like to do instead is to:

1. calculate for each timepoints the fraction of voxels in my ROI where the spikiness measure ‘s’ exceeds a certain threshold. (I think the default threshold in 3dDespike is 2.5)

2. Count for each subject the number of timepoints where this fraction exceeds a certain threshold

Does this sound right and what AFNI function can I use that would allow me to do step 1 and 2?

Thank you very much in advance!
Carolin

[ptaylor]

Hi, Carolin-

Just to zoom out for a second:
A) Are you wanting to characterize the spikiness in an ROI at the start of processing, before anything has been done to the time series?
B) Also, do you want your spikiness measure to reflect temporal associations of spikes? That is, say there are 10 voxels in the left-hippocampus and you have 100 time points. Do you care if 3dDespike finds 10 spikes scattered around at different voxels at different time points, or do you want to know if those 10 spikes occur just within a small number of voxels (so, spatially associated), or do you care if those 10 spikes occurred all in a single time point (temporal association)?

• my guess is that from a QC point of view, having temporal synchronization is the main feature, but that is just a guess—that might be the sign of a motion event or dropout or something, which would suggest that time point’s data might be compromised.
  C) Note that very noisy/problematic time series could have very low spikiness because the variability is so high—so just using spikiness as a measure might mask the problem. So actually, a low variability time series might be susceptible to finding some high spikiness, not because of badness but just because baseline variability is so low. That is why we often use cross-space and within-time checks for badness with outlier fraction count at a given time point—that seems to me to be the most useful role for spike-checking.
  D) How do you plan to use an omnibus spikiness measure to QC your hippocampal data? Do you want to possibly censor out bad voxels, or do you want to possibly censor out bad time points of data, or do you want to know whether a subject can be used at all?

–pt

[carolin31]

Hello Paul,

Sorry for the late reply. I was on a longer break.

Regarding your questions:

A) Yes. I would like to evaluate the spikiness in the HIP before any data processing
B) I guess I am not sure. My thoughts are that I would like to have a QC measure that reflects how ‘bad’ the HIP as an ROI is affected by motion artifacts. From what I understand motion can affect different locations of the brain to a different degree (e.g worse in PFC than in the center of the brain). For my analysis I can not censor any volumes because it requires a temporal continuity of the scans. I am also only interested in the HIP for now. So I would hate to exclude subjects that based on a whole brain motion measure like norm or srms have to be excluded when maybe the motion in the HIP itself is not as bad and the signal in the HIP could be still used for my analysis. Does that make sense to you? Now when you are saying that a temporal synchronization indicates most likely a motion artifact then this is probably something I would want to be considered in my summary measure.
C) I had no idea about the relationship you describe between a low baseline variability in the time series and high spikiness. But it does make a lot of sense to me. In this context I understand better why we typically use the outlier fraction count.
D) As mentioned in B, I can not censor my time series because I need temporal continuity for my analysis. So I would like to use this spikiness measure only as an indicator of how bad the motion artifacts are in the HIP and then to exclude entire subjects when the motion is too bad.

Best
Carolin

[carolin31]

Hello afni experts,

I would like to follow up on my question here in case my post was overlooked.

I still believe that using 3dToutcount is not the right thing to do in this case because outliers are determined based on calculating the trend of ‘s’ across all time points and calculating the trend of ‘s’ does not make sense to me.

To consider the temporal association, which from what I understood in point B would be the best indicator of a motion/signal dropout event, it would make sense to me to calculate for each timepoint the fraction of voxels where s exceeds a certain threshold (e.g. 2.5). I am not sure which afni function I could use for that considering 3dBrickstat would not do this for each timepoint.

Thank you very much in advance
Carolin

[ptaylor]

Hi, Carolin-

I would suspect that removing the trend would actually be useful, because you are looking to check for outliers before any processing; there are baseline fluctuations (low/slow moving trends in the data) that exist, and getting rid of those seems necessary before looking for outliers. Also note that with 3dToutcount, you can choose whether to remove trends or not, by selecting the “-polort …”. The default 3dToutcount command that afni_proc.py creates to count outliers in the brainmask (followed by the 1deval command to turn the count into a censoring list if frac > 5%) is:

3dToutcount -automask -fraction -polort 3 -legendre                     \
    DSET_EPI > FILE_TOUT.1D

1deval -a FILE_TOUT.1D \
     -expr "1-step(a-0.05)" > FILE_CENSOR.1D

To me, it would probably make sense to leave in the “-polort …” option, but that could be removed. I am not sure if you are planning to binarize the output, which is what the 1deval command does here—anything more than 5% turns into a time point to censor.

Because of the issue that outlier checking is inherently related to the scale of “typical” fluctuations in the time series, it is hard not to see that also having a “quantitative” measure of those fluctuations would also be useful. For example, afni_proc.py has the “scale” block whereby each voxel’s fluctuations are mapped with its own baseline and std dev to units of “BOLD % signal change”. You could borrow that idea up front for your unprocessed data, as well. Then, you would see what the typical fluctuations are in your data in terms of BOLD % signal change. So, if you have a couple outliers but the typical fluctuations are miniscule, that would be different than having no outliers but a 50% BOLD signal change on average, say. The code for that looks like (with a mask, optionally):

# ================================= scale ==================================
# scale each voxel time series to have a mean of 100
# (be sure no negatives creep in)
# (subject to a range of [0,200])
3dTstat -prefix DSET_MEAN DSET_EPI
3dcalc -a DSET_EPI -b DSET_MEAN \
           -c DSET_MASK                            \
           -expr 'c * min(200, a/b*100)*step(a)*step(b)'       \
           -prefix DSET_EPI_SCALED

And before estimating local standard deviation or variance, you might want to detrend DSET_EPI_SCALED, say.

3dReHo is a program that calculated Kendall’s Coefficient of Concordance (KCC), which has been dubbed “regional homogeneity” or ReHo in the FMRI world. Conceptually, it states how similar a set of time series are, based on their relative fluctuations, where the set can be >2 (and typically it is a local neighborhood, but can also be an ROI). That might be useful in this question, though there will be degeneracies of interpretations for “high” or “low” values. But it would add another dimension to your investigation.

–pt

[carolin31]

Hello Paul,

Thank you very much for your reply.

I apologize if I repeat my question but I want to make sure that I made clear that I am trying to summarize the amount of spikes in the HIP based on the spikiness measure ‘s’ (output from 3dDespike). I am wondering if using 3dToutcount is appropriate to do so in this case.
So I am not looking at the raw signal of the time series and try to determine outliers but rather I am looking at the spikiness measure ‘s’ from 3dDespike that was calculated by already fitting a smooth-ish curve to the raw signal in each voxel and determining the degree of deviation of each time point under consideration of MAD and standard deviation.
I am not sure conceptually if it makes sense to determine outliers (in this case spikes) in this signal based on calculating again a trend over this spikiness measure and the deviations from this trend. I thought that an outlier (in this case spike) is determined based on the value of s (default in afni >2.5) so it would make more sense to me to now just count the fraction of voxels at each time point where s exceeds lets say 2.5.
Or are you saying that I should detrend the raw time series before running 3dDespike?
Or that what I am trying to do is not a good idea and I should just use 3dToutcount on the raw signal with the -mask option and limit the analysis to the HIP?

I love the idea of also looking and comparing the BOLD % signal change in the HIP voxels as an additional measure for QC. And looking at ReHo. I have not heard of this measure before. Thank you very much for these suggestions!

Carolin

[carolin31]

Dear afni experts,

I wanted to follow-up on my question.
Thank you very much in advance!

Carolin

[rickr]

Hi Carolin,

I am pretty sure Paul had meant running 3dToutcount on the original time series, not the 3dDespike -ssave output. But let me just go from what you have said you want for now. You wanted to:

• run 3dDespike -ssave
• calculate the fraction of voxels in your ROI where s exceeds 2.5, say
• count the number of time points where the fraction exceeds a a threshold

So assuming you have a time series called ssave+orig and a mask mask+orig:

# calculate a time series "mask" of thresholded values
set thr_spike = 2.5
3dcalc -a ssave+orig -expr "step(a-$thr_spike)" -prefix ssave_thresh

# calculate the ROI average at each time point into a 1D time series
3dmaskave -q -mask mask+orig ssave_thresh+orig > ssave_thresh_ave.1D

# then threshold this 1D time series at some max outlier fraction
set thr_ofrac = 0.1
1deval -a ssave_thresh_ave.1D -expr "step(a-$thr_ofrac)" > ssave_outfrac.1D

# then count the number of time points (could grep and wc, or...)
3dTstat -sum -prefix - ssave_outfrac.1D\' > ssave_nt_bad.1D

echo "the number of bad time points is `cat ssave_nt_bad.1D`"

Is this what you are picturing?

• rick

[carolin31]

Hello Rick!

Yes! This is what I envisioned! Thank you very much for your help

Cheers
Carolin

[carolin31]

Hello Rick,

sorry to bother you again with this. As I am working through the code I realize that
step 2 should involve calculating the fraction of spikes (spikes above 2.5 threshold) within the ROI for each time point and not the average.

And step 3 and 4 then to apply a threshold for the fraction and count how many time points exceed this threshold.

What afni function can I use to calculate the fraction within a ROI?

Thank you
Carolin

[rickr]

Since the time series was thresholded into a binary spike time series in step 1, the ROI average is indeed the fraction of voxels that are spikes at each time point. It is not the average spikiness, but the average of the spiked voxels.

For example, if at some time point 40% of the ROI voxels are 1.5 and 60% are 3, the binarized time series will have 40% of the voxels at 0 and 60% at 1, since they are thresholded at 2.5. And then the ROI average will be 0.6, the fraction of voxels that seem “bad”. Note that the values in this average time series are always in [0,1].

Does this seem reasonable?

• rick

[carolin31]

Oh I see! Yes that does make a lot of sense now!

Thank you for your help Rick! :))

Cheers
Carolin

[carolin31]

Hello afni experts,

I have a quick question regarding this topic.

I really liked the idea to also evaluate BOLD % signal change in my data as another quality check for motion.
I scaled the data using the code above and before calculating the stdev I would like to detrend it as recommended. For detrending I used the default setting in 3dTproject (so 3 polynomials). When looking at the detrended signal it looks ok for most voxels but in some I can still see a trend (highlighted in red in the images). Do you recommend to choose a higher polynomial or is there a good reason not too. Again I am processing the signal just to evaluate potential motion artifacts as QC. If I can choose a higher number for polynomials how do I know how high this number can be? Do I just keep running the analysis and checking the signal until I am satisfied with the detrended signal?

Thank you very much
Carolin