Open julfou81 opened 3 years ago
Hi again,
I went on trying to understand what is going on and I believe that there is a weird interaction between the SBref volume taken as the reference for the realignment for the bold run and the cost function used in mcflirt. For information, the TR we used is 1.1s and the contrast in the SBref image is much stronger that for the volumes in the bold run. This is expected and that is why the SBref volume is used to realign the bold space to the anatomical space usually with some boundary based registration. But for mcflirt, there may be an issue and the difference of contrast between the SBref volume and a bold volume is not the same as the difference of contrast between two bold volumes from the same run. Here is the SBref volume: and the first volume of the corresponding run:
To illustrate my point, here is an observation: the plot of the second column of the .par file, output by mcflirt, with three different execution related to the same bold run:
=> there is something fishy with mcflirt, SBref and the cost function normcorr ! I never noticed that before even though we used fmriprep for many projects and subjects, maybe it is especially problematic with subjects with extremely low movements? I will also ask the FSL forum about this point as it is an issue not specific to fmriprep but rather to FSL-mcflirt. Also since this approach SBref - mcflirt is also used by the HCP, I will dig in that discussion to see if that was observed already.
To summarize: the mcflirt behavior (giving a lot of 0 in the motion parameters) seems to be coming from three factors:
Anyhow, even though the motion corrected runs seem fine, I think it is more problematic to have a lot go zeros, than to have none. In some runs, one or two columns of the motion parameters are only zeros, and therefore so are the derivative and squared ones, which is problematic if one uses those columns full of zeros in the GLM as nuisance regressors. This affects also the Framewise displacement estimation, one can see the difference between the FD curve calculate by MRIQC (with AFNI 3dvolreg) and the one calculated with with FMRIPREP with FSL- mcflirt:
FMRIPREP (mcflirt)
MRIQC(3dvolreg)(same run as above)
MRIQC (mcflirt) (same run)
FD_mean calculated by FMRIPREP: 0.035mm FD_mean calculated by MRIQC (AFNI 3dvolreg): 0.112 mm FD_mean calculated by MRIQC (FSL mcflirt): 0.104 mm
bids_name | fd_mean_mriqc_3dvolreg | fd_mean_mriqc_mcflirt | fd_mean_fmriprep |
---|---|---|---|
sub-testpilote2_ses-01_task-VibStim_run-01_bold | 0,194680121 | 0,178982256 | 0,136 |
sub-testpilote2_ses-01_task-VibStim_run-02_bold | 0,186500352 | 0,175467697 | 0,136 |
sub-testpilote2_ses-01_task-VibStim_run-03_bold | 0,196197217 | 0,180603697 | 0,146 |
sub-testpilote2_ses-01_task-VibStim_run-04_bold | 0,18933964 | 0,174164599 | 0,143 |
sub-testpilote3_ses-01_task-VibStim_run-01_bold | 0,116184151 | 0,100634748 | 0,06 |
sub-testpilote3_ses-01_task-VibStim_run-02_bold | 0,103993465 | 0,090084196 | 0,064 |
sub-testpilote3_ses-01_task-VibStim_run-03_bold | 0,097314613 | 0,092405695 | 0,038 |
sub-testpilote3_ses-01_task-VibStim_run-04_bold | 0,105072874 | 0,111065972 | 0,013 |
sub-testpilote3_ses-02_task-VibStim_run-01_bold | 0,112095188 | 0,104452559 | 0,035 |
sub-testpilote3_ses-02_task-VibStim_run-02_bold | 0,112124946 | 0,11514835 | 0,025 |
sub-testpilote3_ses-02_task-VibStim_run-03_bold | 0,115953908 | 0,114734764 | 0,021 |
sub-testpilote3_ses-02_task-VibStim_run-04_bold | 0,114908759 | 0,112333997 | 0,015 |
sub-testpilote3_ses-02_task-VibStim_run-05_bold | 0,114918923 | 0,102249177 | 0,018 |
I went to plan B for this case of low movers where obviously the motion estimation with the combo mcflirt - SBref - cost function (normcorr) was wrong (too many zeros). For plan B, I decided to not use the SBref images anymore, so FMRIPREP would have to build its own reference scan by averaging non-steady state volumes for each run. I am happier with the results this way. No more zeros in the motions parameters and FD value are much more similar to the results from 3Dvolreg or mcflirt implemented by MRIQC:
bids_name | fd_mean_mriqc_3dvolreg | fd_mean_mriqc_mcflirt | fd_mean_fmriprep | fd_mean_fmriprep_noSBref |
---|---|---|---|---|
sub-testpilote2_ses-01_task-VibStim_run-01_bold | 0,194680121 | 0,178982256 | 0,136 | 0,175 |
sub-testpilote2_ses-01_task-VibStim_run-02_bold | 0,186500352 | 0,175467697 | 0,136 | 0,164 |
sub-testpilote2_ses-01_task-VibStim_run-03_bold | 0,196197217 | 0,180603697 | 0,146 | 0,174 |
sub-testpilote2_ses-01_task-VibStim_run-04_bold | 0,18933964 | 0,174164599 | 0,143 | 0,165 |
sub-testpilote3_ses-01_task-VibStim_run-01_bold | 0,116184151 | 0,100634748 | 0,06 | 0,095 |
sub-testpilote3_ses-01_task-VibStim_run-02_bold | 0,103993465 | 0,090084196 | 0,064 | 0,09 |
sub-testpilote3_ses-01_task-VibStim_run-03_bold | 0,097314613 | 0,092405695 | 0,038 | 0,08 |
sub-testpilote3_ses-01_task-VibStim_run-04_bold | 0,105072874 | 0,111065972 | 0,013 | 0,082 |
sub-testpilote3_ses-02_task-VibStim_run-01_bold | 0,112095188 | 0,104452559 | 0,035 | 0,096 |
sub-testpilote3_ses-02_task-VibStim_run-02_bold | 0,112124946 | 0,11514835 | 0,025 | 0,087 |
sub-testpilote3_ses-02_task-VibStim_run-03_bold | 0,115953908 | 0,114734764 | 0,021 | 0,086 |
sub-testpilote3_ses-02_task-VibStim_run-04_bold | 0,114908759 | 0,112333997 | 0,015 | 0,082 |
sub-testpilote3_ses-02_task-VibStim_run-05_bold | 0,114918923 | 0,102249177 | 0,018 | 0,084 |
There must be more validation for this case but eventually, if this behavior of mcflirt is confirmed on low movers, what do you think of a pull request to create a new heuristic: when subjects have low motion and mcflirt results are too low with SBref as reference, then drop the SBref image and move to the next method to calculate a reference image? Or give an option to theuser to change the cost function in mcflirt to "mutualinfo" which seens to do the job in my case, while keeping the SBref image as a target.
As this strange behavior during HMC seems to be related to mcflirt, data were sent to the FSL mailing list. FSL developers reproduced the bug internally, they are investigating. Meanwhile looking at #2367, I wonder if the behavior we observe (different HMC results depending on using SBref as reference image or not) is somehow linked to a difference in the qform and sform between SBref and the related MB images. I will check that also. For us the images are acquired using the C2P from CMRR which is widely used, and converted using dcm2niix, I am doubtful that there is something there but I will check. As an illustration, here is excerpt from the FMRIPREP report for the same run, whether the SBref images were present (on top) in the func folder or if they were removed (on bottom). Notice how the FD trace change from one condition to the other. FD mean goes from 0.05mm (SBref) to 0.09mm (no SBref). MRIQC for the same run (with AFNI 3dvolreg) was giving FDmean = 0.11mm.
EDIT: qform and sform matrices are fine in my case. Orientation information is totally identical between SBref images and corresponding bold run.
EDIT: FSL team is well aware of this problem, they were able to reproduce the mcflirt behavior with sbref images as target file for motion correction but no solution is coming so far. The HCP team also encountered this problem a few years ago and it is something that was not solved at the time either. When SBREF images are involved, could we think of changing the workflow and do something similar to what is done with the HCPpipelines? i.e using the SBref image to compute the transformation from bold space to T1w space, use the mean of the 10 volumes (or mean of non steady state volumes) for the motion correction target and compute the transformation from this mean to the SBref image. And concatenate all the the transformations: motion correction matrix for each volume - transfo from mean to SBref - transfo from SBref to T1w. I am new with nipype but I could help on this PR if one think that it would be a good direction to take.
Must be related to : #1684
I just checked: the problem with mcflirt is not only related to a mismatch between orientation header between the SBref image and the related multiband serie , as it is mentioned in #1684.
Example: Two files:
Difference from the two headers ( diff
command between the two outputs of fslhd
)
1c1
< filename sub-pilot2_task-PhantomMobile_run-05_sbref.nii.gz
---
> filename sub-pilot2_task-PhantomMobile_run-05_bold.nii.gz
5c5
< dim0 3
---
> dim0 4
9c9
< dim4 1
---
> dim4 387
34,35c34,35
< slice_name sequential_decreasing
< slice_code 2
---
> slice_name Unknown
> slice_code 0
66,67c66,67
< descrip TE=35;Time=114425.535;phase=1
< aux_file Single-band_reference_S
---
> descrip TE=35;Time=114430.313;phase=1;mb=4
> aux_file Unaliased_MB4_PE3_LB_SE
=> The orientations are strictly identical. However here is the FD trace with FMRIPREP with or without --ignore sbref
Same run, preprocessed with fmriprep v20.2 with the SBref file use as reference for mcflirt and bold-to-T1w registration according to the fmriprep heuristic on choice for the bold reference:
Same run, this time preprocessed with fmriprep v20.2.2 with the SBref discarded by the --ignore sbref
option:
FD mean goes from 0.05 mm to 0.13 mm!
Hi, apologies for not getting back to you sooner. I've had a skim of what you have above, and appreciate your thorough reports. I'm honestly not sure what to make of this, but the notion that it's an issue with MCFLIRT's cost function seems likely to be correct.
Just to clarify (and apologies if I missed this) have you visually compared the motion-corrected time series, e.g. by animating the mid-sagittal plane) to verify if both are equally well-corrected, or if one works more poorly?
One additional thing to try is to see whether AFNI's 3dvolreg has similar issues or is more robust to the choice of reference. In the long run, we want to move to that from MCFLIRT if we can.
Hi, thank you for your message. Don't apologize, I raised that issue to check if somebody else saw this, and to see if it can be improved.
I carried out the comparison you asked: visually, looping through the two corrected time series in fsleyes (from fmriprep with --ignore sbref
option or not), both look very similar.
I did carried out the extra test with AFN's 3dvolreg just now, either by taking the SBref as reference, or the first volume, and looking at the motion parameters, I don't see any zeros and they look pretty similar.
Again looping through all these time series (I re-ran MCFLIRT alone as well, to have a fair comparison), I have a visual impression that 3dvolreg (that I never used before this test) gives better motion correction that MCFLIRT, either with SBref or the first volume as target for the realignment.
I will provide more quantitative measurements later from this comparison.
This is one single test, one subject, one run, there must be of course more validation, but at first glance, I would encourage you to follow your idea to go with 3dvolreg for motion correction in FMRIPREP.
Regarding the mcflirt - 3dvolreg comparison (outside the problem with SBref), both are available in MRIQC, did you get any return showing that one was more robust the other?
Here is one plot which is illustrative: These are the 3 translation parameters for the same run, for which the target for realignment was the SBref image, with either 3dvolreg or mcflirt:
Hello experts, thank you for this intriguing discussion. We use the data from HCP-Aging cohort and observe that the motion regressors from some subjects (say "some" is very conservative actually... we ran 12 subjects so far and more than half of the 12 subjects have at least one resting functional run that shows zeros in one or more motion directions). I found that this issue is also reported on the HCPpipeline github discussion (thank you @julfou81 !!).
From my understanding, the alternatives include (with the sbref as reference) (1) use "mutual information" cost function (FSL mcflirt) (2) use AFNI 3dvolreg: this option seems to be more recommended by the above discussion
We are currently at the stage of deciding whether to use the HCP-Aging minimally preprocessed data or to run the preprocessing on our own with fMRIPrep, therefore, we are comparing the outputs from two pipelines. The zeros in motion regressors appear in low-motion runs from both HCP minimally preprocessed regressor and fMRIPrep confound output.
I am wondering what is the approach that you would recommend to proceed with this issue?
Any opinions or suggestions will be very helpful! Thank you.
Thank you @yenwchen for reporting also about this issue.
This issue is in general very visible for very "good" subjects, i.e. subjects that do not move very much during the FMRI runs.
Here is a new case that I came across, where the difference in motion estimation is rather quite different (especially at the beginning of the run) if I compare MRIQC (using 3dvolreg), FMRIPREP v21.0.2 (with SBref used as reference) or FMRIPREP v21.0.2 (with --ignore sbref
):
MRIQC:
FMRIPREP v21.0.2 with SBref used as reference:
FMRIPREP v21.0.2 with --ignore sbref
:
I would suggest to move toward adding an option in FMRIPREP to automatically change the cost fonction of mcflirt
when SBref is used, which is a relative quick fix before moving to 3dvolreg
.
I didn't read all of that...but years ago I was writing my own confound and detrend program. I found that the -plots
output from mcflirt
didn't make sense to me. I used the matrix outputs from -mats
and decomposed them using the Euler angle module from Nibabel. Worked great.
What version of fMRIPrep are you using?
FMRIPREP v 20.2.1
What kind of installation are you using? Containers (Singularity, Docker), or "bare-metal"?
Singularity installation (version 2.6.1-dist), running on a high performance computing cluster.
What is the exact command-line you used?
Have you checked that your inputs are BIDS valid?
yes the input passes the bids-validator
Did fMRIPrep generate the visual report for this particular subject? If yes, could you share it?
Yes the visual report was generated, no error reported, it is under the dropbox link at the end of this post.
Can you find some traces of the error reported in the visual report (at the bottom) or in crashfiles?
no error reported
Are you reusing previously computed results (e.g., FreeSurfer, Anatomical derivatives, work directory of previous run)?
No
fMRIPrep log
If you have access to the output logged by fMRIPrep, please make sure to attach it as a text file to this issue.
We observed this strange output while building the GLM with nuisance regressors and noting that the columns with zeros where correlated with each other.
output: HTML reports, output and confounds_timeseries.tsv files under the following link: https://www.dropbox.com/sh/3e455qbh8anqwl3/AAACn2-tnmRBhXx44M_vO0qEa?dl=0
Here is for example the content of the .par file output of mcflirt run by fmriprep: the zeros are very strange and it is a very different result from what I obtain if I use SPM for realigning the same run (but realignment is to the first scan instead to the sbref as it is set with fmriprep) No zeros are found with SPM. EDIT: I just check with my local FSL installation and I found that mcflirt is giving the same output as fmriprep! Thus it is not seem to be a problem with fmriprep, rather a strange behaviour of mcflirt! This subject had very low movements. Thank you for sharing your thoughts!
Here is a plot (made with FSLeyes) of the translation along x, calculated with mcflirt (.par file, 4th column, ref= sbref) and with spm (rp file, 1st column, ref=1st scan): Quite different! Another example: same run, looking at the ortation around y, calculated with mcflirt (.par file, 2nd column, ref= sbref) and with spm (rp file, 5th column, ref=1st scan):