Question re: filtering FD values

[Danella1]

Hi again, xcp-d experts -

Thanks for all of your help thus far! I now have a question about filtering out respiratory effects from FD values. I have been influenced by papers from Damien Fair’s and Caterina Gratton’s labs to try this. I have tried both notch and lp filters in xcp-d runs - and I am expecting to see the FD traces in the QC plots decrease, as they do in the papers mentioned above (comparing .html files from fmriprep and xcp-d). However, the FD traces that I see in these outputs are identical. Also, this filtering is not mentioned in the xcp-d methods boilerplate, so I am wondering if the filtering is actually taking place? How can I check? And if it indeed ran, would there be any way to have the filtered FD values in the output .html, to be able to see the effect of this filtering?

Here is the code that I have run:

singularity run --cleanenv -B /ix/dhafeman:/home/xcp_d -B /ix/dhafeman/license.txt:/opt/freesurfer/license.txt /ix/dhafeman/singularity/xcp_d-0.1.3.simg /home/xcp_d/BDLONG/Preproc/fmriprep /home/xcp_d/BDLONG/Preproc participant -f 0.1 -w /home/xcp_d/BDLONG/Work --smoothing 8 --despike -p 36P --lower-bpf .009 --upper-bpf .08 --motion-filter-type lp --participant_label 400

and

singularity run --cleanenv -B /ix/dhafeman:/home/xcp_d  -B /ix/dhafeman/license.txt:/opt/freesurfer/license.txt /ix/dhafeman/singularity/xcp_d-0.1.3.simg /home/xcp_d/BDLONG/Preproc/fmriprep /home/xcp_d/BDLONG/Preproc participant -f 0.1 -w /home/xcp_d/BDLONG/Work --smoothing 8 --despike -p 36P --lower-bpf .009 --upper-bpf .08 --motion-filter-type notch --band-stop-min 12 --band-stop-max 20 --participant_label 400

[Steven]

Hi Danella,

@Danella1

I have tried both notch and lp filters in xcp-d runs - and I am expecting to see the FD traces in the QC plots decrease

This doesn’t seem intuitive to me, as FD should be more of a function of the inter-volume registration (that is, about head geometry and gross head movement), while the filters are applied to the signal, irrespective of geometry. Thus, I would expect the FD traces to be the same, but the fMRI signal to be different for filter vs. no filter.

Best,
Steven

[tsalo]

@Steven, I believe she’s referring to the motion parameter filtering, rather than the BOLD bandpass filtering.

@Danella1 I think there might be a bit of a mixup with the parameters. I think you want to set --band-stop-min and --band-stop-max. Otherwise, the motion filtering step may not actually filter the motion parameters, even though you selected the “lp” filter type. Please note that these parameters are in breaths-per-minute, rather than Hertz. There are some recommended values here, based on participant age.

singularity run --cleanenv \
    -B /ix/dhafeman:/home/xcp_d \
    -B /ix/dhafeman/license.txt:/opt/freesurfer/license.txt \
    /ix/dhafeman/singularity/xcp_d-0.1.3.simg \
    /home/xcp_d/BDLONG/Preproc/fmriprep \
    /home/xcp_d/BDLONG/Preproc \
    participant -f 0.1 -w /home/xcp_d/BDLONG/Work \
    --smoothing 8 --despike -p 36P \
    --lower-bpf .009 --upper-bpf .08 \  # these determine the BOLD bandpass filter
    --motion-filter-type lp \  # this enables the motion filter, but the default filter values are 0 and 0 (i.e., it is skipped)
    --participant_label 400

It’s also possible that the HTML report simply loads the wrong time series when plotting FD. You should find the filtered FD in the *_desc-framewisedisplacement_bold.tsv file (renamed to *_desc-framewisedisplacement_motion.tsv in the upcoming release).

[Danella1]

Thanks for the information. Indeed, the *_desc-framewisedisplacement_bold.tsv shows lower values than those from fmriprep. I am wondering about the lp defaults though. I saw on github that the default for the lp was .1 Hz (xcp_d/confounds.py at 33e6fa843d39db27da571c5f1a6b4396a25c57d3 · PennLINC/xcp_d · GitHub). Is that incorrect?

If I were to set it manually (using lp), would I set --band-stop-max 6?

Thank you for your help with this!

[tsalo]

@Danella1

I saw on github that the default for the lp was .1 Hz

That is a completely reasonable interpretation. Unfortunately, the documentation is incorrect (probably my fault). cutoff is actually in breaths per minute, and that default value is never used, because it’s taken directly from the band-pass values. Within that function, the cutoff and band-pass values are all divided by 60 to convert the breaths-per-minute to Hertz before applying the filter.

@Danella1

If I were to set it manually (using lp), would I set --band-stop-max 6?

I believe that should work. However, 6 might be too low for filtering the parameters, based on the documentation.

[Danella1]

Ok I will rerun it like that. Gratton’s group has found <0.1 Hz to work well, so I wanted to start there (Removal of high frequency contamination from motion estimates in single-band fMRI saves data without biasing functional connectivity - ScienceDirect). I realize this is different from the notch filtering parameters - and to be honest, I have a fairly superficial grasp of this, so please let me know if you think I am misunderstanding. Thank you!

[tsalo]

Sorry, the confusion was entirely on my end. I completely forgot that the “notch” option is a band-stop filter instead of a band-pass filter. 0.1 Hz/6 BPM makes sense. It also means that you would preferably use --band-stop-min instead of --band-stop-max to set this value, but in the current release (0.1.3), that doesn’t matter since it uses the minimum value between --band-stop-min and --band-stop-max. In future releases, there will be an error if folks use --motion-filter-type lp with --band-stop-max.

Actually… I worry that setting --band-stop-max without setting --band-stop-min might resolve to zero when you run it, since the low-pass filter’s cutoff will be determined from numpy.min([band_stop_min, band_stop_max]), and the one you don’t set will default to 0… Please follow up if that seems to be the problem.

[Danella1]

Thank you! BTW, I ran lp with and without --band-stop-min 6, and got the same answer (which was similar to the notch filter, and substantially less than the unfiltered FD). So, I think this might be the default value? Thanks!

Two additional questions: Any chance that the .html output can reflect the filtered FDs in a future iteration? Also, is there anywhere in the xcp_d output where it is described whether/how FDs were filtered? It would be helpful if that were integrated in the methods boilerplate, or if there is some other way to know from the results.

Thanks again for all your help with this!

[tsalo]

@Danella1

I ran lp with and without --band-stop-min 6, and got the same answer (which was similar to the notch filter, and substantially less than the unfiltered FD). So, I think this might be the default value?

Oh wow, that’s an annoying bug. It looks like v0.1.3 completely ignores the --band-stop-min and --band-stop-max parameters when the motion-filter-type is set to ‘lp’, and just uses 0.1 no matter what. Unfortunately, that 0.1 is still divided by 60 as if it was in BPM, rather than Hertz, so the low-pass filter in v0.1.3 appears to be unusable.

When I said that it uses the minimum value between the two parameters, I was looking at more recent code, and thus my advice was completely wrong. I’m really sorry about this!

@Danella1

Any chance that the .html output can reflect the filtered FDs in a future iteration?

Yes, absolutely. It looks like the plots actually all show the filtered FD rather than the original, but I think the key is being able to see both before and after filtering. I’ve opened an issue about this: Preprocessing QC plot only shows FD after filtering · Issue #576 · PennLINC/xcp_d · GitHub.

@Danella1

Also, is there anywhere in the xcp_d output where it is described whether/how FDs were filtered? It would be helpful if that were integrated in the methods boilerplate, or if there is some other way to know from the results.

The filtering isn’t reflected in the boilerplate at the moment. I have opened a pull request to add this information, and we’ll make sure to merge it before our next release, which will hopefully be soon.

[Danella1]

Thank you!

I am a bit confused though. It really does seem that the lp_default is doing something (since it is very similar to notch, and different from fmriprep/no_filter. Do you think a filter at .1/60=.001667 could have this effect?

Also, here are my plots for fmriprep and xcp (after filtering - this is with notch). They seem identical to me? I am not sure why this would not be showing the filtered FDs:
fmriprep:

xcp_d:

I really appreciate all of your help with this!

[tsalo]

They definitely do look identical in the line plots, but the values you show in the Excel file are different. My guess is that the xcpd report is showing the wrong time series.

However, I have opened a pull request to fix this: Plot appropriate version of FD in pre- and post-processing figures by tsalo · Pull Request #579 · PennLINC/xcp_d · GitHub. Here is what the updated pre- and post-processing plots look like when using a notch filter from 12 to 18 BPM on test data:

As you can see, in the new version, the filtered time series will be plotted in the post-processing figure, while the unfiltered one will be shown in the pre-processing figure. Hopefully we’ll be able to release the new version by the middle of next week.

[Danella1]

Wonderful! I will wait to run my analyses until the new release. Thank you, and have a good weekend!

[tsalo]

Hi @Danella1,

It was correctly pointed out to me that, since the filtered framewise displacement is just used to flag outlier volumes for censoring, and not for the nuisance regression, it would be misleading to use the filtered framewise displacement in the postprocessing carpet plot. Instead, I’ve added a censoring plot to the report, like the one below.

Also, we’ve released 0.2.0, so I’m hoping it will work better for you.

[Danella1]

Thank you! The output and reports in version 0.2.0 look great!

[Danella1]

I spoke too soon. xcp_d worked very well for a single task, but when I ran my first subject, I had several crash reports. They were all essentially the same, and only related to the node xcpd_wf.single_subject_400_wf.nifti_postprocess_2_wf.execsummary_wf.plot_svgx_wf. It seems that everything else ran just fine - there is just no Executive Summary .html. There is not really a pattern to which tasks led to the error - they were pretty evenly distributed across sessions, and different tasks for each session (e.g. ses-1: movie, rest-run2; ses-2: inscapes_run2, rest_run2, reward).

EDITED TO ADD: I figured out what is going on. The exec summary node runs fine when there is no censoring. However, when there is censoring, the “regressed_data” input does not match the length of “raw_data” or “filtered_motion”, and thus gives the error. Please let me know if I should post on github! Thank you!

Additional info below…

Here is the crash message. I would like to sort this out before running all of my subjects, because I do find the Executive Summary quite helpful.

Node: xcpd_wf.single_subject_400_wf.nifti_postprocess_2_wf.execsummary_wf.plot_svgx_wf
Working directory: /home/xcp_d/BDLONG/Work/xcpd_wf/single_subject_400_wf/nifti_postprocess_2_wf/execsummary_wf/plot_svgx_wf

Node inputs:

TR = 1.5
filtered_motion =
mask = /home/xcp_d/BDLONG/Preproc/fmriprep/sub-400/ses-01/func/sub-400_ses-01_task-movie_space-MNI152NLin6Asym_desc-brain_mask.nii.gz
rawdata = /home/xcp_d/BDLONG/Preproc/fmriprep/sub-400/ses-01/func/sub-400_ses-01_task-movie_space-MNI152NLin6Asym_desc-preproc_bold.nii.gz
regressed_data =
residual_data =
seg_data =

Traceback (most recent call last):
File “/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/plugins/multiproc.py”, line 67, in run_node
result[“result”] = node.run(updatehash=updatehash)
File “/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py”, line 527, in run
result = self._run_interface(execute=True)
File “/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py”, line 645, in _run_interface
return self._run_command(execute)
File “/usr/local/miniconda/lib/python3.8/site-packages/nipype/pipeline/engine/nodes.py”, line 771, in _run_command
raise NodeExecutionError(msg)
nipype.pipeline.engine.nodes.NodeExecutionError: Exception raised while executing Node plot_svgx_wf.

Traceback:
Traceback (most recent call last):
File “/usr/local/miniconda/lib/python3.8/site-packages/nipype/interfaces/base/core.py”, line 398, in run
runtime = self._run_interface(runtime)
File “/usr/local/miniconda/lib/python3.8/site-packages/xcp_d/interfaces/surfplotting.py”, line 206, in _run_interface
‘after_process’] = plot_svgx(
File “/usr/local/miniconda/lib/python3.8/site-packages/xcp_d/utils/plot.py”, line 523, in plot_svgx
DVARS_timeseries = pd.DataFrame({
File “/usr/local/miniconda/lib/python3.8/site-packages/pandas/core/frame.py”, line 529, in init
mgr = init_dict(data, index, columns, dtype=dtype)
File “/usr/local/miniconda/lib/python3.8/site-packages/pandas/core/internals/construction.py”, line 287, in init_dict
return arrays_to_mgr(arrays, data_names, index, columns, dtype=dtype)
File “/usr/local/miniconda/lib/python3.8/site-packages/pandas/core/internals/construction.py”, line 80, in arrays_to_mgr
index = extract_index(arrays)
File “/usr/local/miniconda/lib/python3.8/site-packages/pandas/core/internals/construction.py”, line 401, in extract_index
raise ValueError(“arrays must all be same length”)
ValueError: arrays must all be same length

Also, in case it is helpful, here was the code that I ran on our computing cluster:

singularity run --cleanenv -B /ix/dhafeman:/home/xcp_d -B /ix/dhafeman/license.txt:/opt/freesurfer/license.txt /ix/dhafeman/singularity/xcp_d-0.2.0.simg /home/xcp_d/BDLONG/Preproc/fmriprep /home/xcp_d/BDLONG/ participant -f 0.1 -w /home/xcp_d/BDLONG/Work --smoothing 8 --despike -p 36P --lower-bpf .009 --upper-bpf .08 --motion-filter-type lp --band-stop-min 6 --participant_label 400

Please also let me know if you would like me to post this as a new issue on neurostars or github. Thank you in advance for your help!

[tsalo]

There is a known bug in that part of the code, but I thought we had a workaround in place. It would be great if you could open an issue on GitHub. Thanks!

[tsalo]

The issue is Crash on execsummary_wf.plot_svgx_wf when there is censoring · Issue #612 · PennLINC/xcp_d · GitHub (for posterity).

[Danella1]

Thank you for fixing this! Will this change be in version “0.2.0” now, or should I use “latest”? Thanks!

[tsalo]

You should use latest for now. The change will be in 0.2.1, which we will hopefully release soon, but not 0.2.0. Also, props go to @kahinimehta for actually solving the problem.

[Danella1]

Hi @kahinimehta and @tsalo, thanks for your help with this! I tried downloading “latest” from singularity, and then reran the analysis - but I am still getting the same error message. Do you have any additional suggestions? I can also paste in the crash message, but it is essentially the same as before. Thank yoU!

It looks like I am running version xcp_d version 0.2.0+0.g46c71f7.dirty.

I also posted this on github (#612), but I wasn’t sure if that message would go through, since the issue is closed.

Thanks!