translationalneuromodeling / tapas

TAPAS - Translational Algorithms for Psychiatry-Advancing Science
https://translationalneuromodeling.github.io/tapas/
GNU General Public License v3.0
213 stars 89 forks source link

Noise ROI Extraction: Erosion, PCA components, and smoothed vs unsmoothed data #66

Closed mrikasper closed 2 years ago

mrikasper commented 5 years ago

From a user:

[...] I am currently analysing my task-related fMRI data.

In my first GLM I only include the 6 movement parameters as regressors of no interest. However, I would like to see whether I can reduce the amount of noise in my statistical maps by including the mean time series for WM and CSF. So I want to include it at the first level (so after the preprocessing). The reason is that some of the noise coming from WM or CSF overlaps or spreads to regions of interest in grey matter (my main focus lies on hippocampus and striatal regions).

So you want to add the WM and CSF time series regressors as additional columns to your GLM design matrix on the first level. Perfect, this is the standard use case of PhysIO.

I have implemented all the steps you described (biggest part was described here: https://www.nitrc.org/forum/forum.php?thread_id=9959&forum_id=2) and I was able to extract the mean timeseries with the PhysIO toolbox.

Great. Points (1)-(4) in the post are generic to both resting state and task-based fMRI. I just wanted to highlight here that the filtering out via writing the residuals described later in that post is not necessary (and even detrimental) for the task-based analysis, because we want to model effects of interest and noise together in order to not mis-attribute shared variance. From what you write, you didn't intend to do that, I just wanted to clarify, because the question comes up frequently.

However there is still something I am not sure about and this is whether I should use the smoothed or unsmoothed normalized functional images. I ran the pipeline on both and the output seems very similar, but I would like to be sure. So far the literature has not really helped me because it's mainly focused on resting-state fMRI.

Would you advise to run it on the smoothed or unsmoothed normalized epi images?

If you, as you mentioned above, are concerned about the proximity of WM/CSF and your regions of interest (Hippocampus, striatum), in that physiological noise could be injected into GM regions, you should also be the other way around, i.e., that valid activation from GM ends up in the WM/CSF "noise regions". Therefore I would advise to use the unsmoothed normalized images for the extraction here. You might lose some SNR on identifying the noise-driven mean time series, but this is not as bad as including signal of interest in the regressor and regressing it out erroneously.

Additionally, the parameter you need to set with respect to cropping, could you elaborate a bit on this (is this related to the eroding of the tissue maps)?

Exactly. Starting from the thresholded tissue probability map for WM/CSF, the erosion removes the amount of voxels at the border of the ROIs that are specified. This is another way of ensuring you end up with WM/CSF voxels only, as opposed to voxels with a partial GM volume.

Lastly, I was not planning to include additional components (PCA), what would you advise with respect to this point?

Theoretically, the PCA components allow to describe fluctuations of similar dynamics between voxels that might have opposing signs (or phases). Imagine you have a perfect sin(t) in one voxel and -sin(t) in another one, and that's your 2-voxel ROI. The mean would just be flat, but there is arguably some shared dynamics in the ROI, and the PCA would determine that, and identify the sin(t) as relevant time series to explain the fluctuations in the ROI. For physiological processes like cardiac pulsatility that have a common driver, but maybe regionally different phase expression, this seems like a useful way to extract the driver.

In practice, try out to add some PCA components and see how much variance they explain (You can also set a component number to sth smaller than 1, e.g., 0.9 to get as many components as needed for removing 90% of the variance). If the first components only explain a small percentage of the total variance, then there is not much structure in the data and you don't have to include them in the model. If they explain, e.g., 30 or 50% of the variance, you can try including them and have a look, if the explain significant variance via an F-test in the 1st level GLM on just that regressor.

My current model includes already 3 regressors of interest, 6 regressors of no interest (movement) and the average timeseries for the 2 noise ROIs per session, and I have 2 sessions .

In terms of model complexity, the number of reasonable regressors without a risk of overfitting depends highly on the number of volumes in your EPI time series. 11 regressors per session does not sound like a lot to me, if you have a couple of hundred volumes.

I hope that helps!

All the best, Lars

NinaDolfen commented 5 years ago

Dear Lars,

thank you for this elaborate reply! Your post answered all my questions. As you say I will explore the PCA options and see whether my model benefits from it or not. I include on average more than 400 volumes per session, so I guess I should not worry about overfitting.

I have one follow up question. In the segmentation step in the preprocessing pipeline I already opt to erode the tissue maps. In the segmentation step it's referred to as 'Warp - Clean up' and you can choose between no clean up, light clean up and thorough clean up. I use the light clean up in accordance with what colleagues in the lab did before. (note that I don't have more information on where this different options stand for, all the information I have on the clean up procedure comes from the manual).

Given the information above, do I need to do additional cropping during the noise ROI extraction procedure? If yes, what is the default value for this parameter?

Thank you in advance,

Nina

NinaDolfen commented 5 years ago

Dear Lars,

Maybe a final follow-up question. I keep getting the following warning: Warning: fMRI volume and noise ROI have different orientation : Afther this the tissue maps are coregistered to the functional.

This is very suprising to me because in my preprocessing I normalize the functional and structural images to my Dartel template in MNI space. Next I segment the normalized structural to get 'normalized' tissue maps. These are used for the extraction of the noise ROI signal. When I visually check the registration between normalized functional and the c2/c3.nii it seems to be fine.

Is this something I should worry about or is it fine as the additional coregistration mainly has positive effects and improves the outcome of the extraction procedure?

Thank you for your input!

Nina

benoitberanger commented 5 years ago

Dear Nina,

In order to project the ROI to the fMRI volume, the ROI volume needs to have the same "orientation" as the fMRI volume. It means same transform matrix (orientation in space), same voxel dimension and same matrix dimension. Both volumes need to have the exact same "grid".

TAPAS checks if the grids are the same, and if not, performs a Coregister : Realign & Reslice. Then it continues with the freshly created noiseROI volume in fMRI volume "grid".

Hope this helps you.

Best, Benoît

NinaDolfen commented 5 years ago

Hi Benoit,

thanks for your reply! I am aware of the aim of the coregistration procedure but given my preprocessing pipeline this should be taken care of before. So the tissue maps are supposed to be coregistered to the functional images.

I extracted the average timeseries and subsequently ran the first level model for 2 different subjects (GLM with 6 movement regressors and 2 noise ROIS). Unfortunately the addition of noise ROI timeseries induces more noise in my maps. I think I must have overlooked something. I played with the threshold but this doesn't seem to influence much (used .9, .95 & .99).

I expected to either see no change in the individual maps at the first level, or in the best case improvement (but not more noise).

Any suggestions?

Nina

benoitberanger commented 5 years ago

Dear Nina,

Just to be sure : 1) Did you use the normalized non-smoothed data ? 2) Did you use intensity threshold and/or cropped voxel to ensure your final ROI is inside the compartment and there is no partial volume ?

Best, Benoît

NinaDolfen commented 5 years ago

Dear Nina,

Just to be sure :

  1. Did you use the normalized non-smoothed data ?

Indeed, 100% sure I used the normalized non smoothed functional images.

  1. Did you use intensity threshold and/or cropped voxel to ensure your final ROI is inside the compartment and there is no partial volume ?

For the first part of your question I am not sure what you are refering to. Do you mean the bias regularisation that takes place during the segmentation? This is done. Or do you refer to the threshold parameter you can set in the toolbox? I might have misinterpreted the use of this parameterr, see below the values I use.

     matlabbatch{nRun}.spm.tools.physio.model.noise_rois.yes.force_coregister = 'No';
    matlabbatch{nRun}.spm.tools.physio.model.noise_rois.yes.thresholds= 0.95;
    matlabbatch{nRun}.spm.tools.physio.model.noise_rois.yes.n_voxel_crop = 0;
    matlabbatch{nRun}.spm.tools.physio.model.noise_rois.yes.n_components = 0;

With respect to the cropping, I set it to 0 for now, see my reply to Lars above where I ask whether the eroding (or clean-up) that takes place in the segmentation step is sufficient or whether I need to do additional cropping (the second message in this issue).

Thank you!

Best, Benoît

benoitberanger commented 5 years ago

I'm referring to PhysIO parameters :

matlabbatch{nRun}.spm.tools.physio.model.noise_rois.yes.thresholds= 0.95;
matlabbatch{nRun}.spm.tools.physio.model.noise_rois.yes.n_voxel_crop = 0;

Step 1 : thresholds = 0.95 A 95% intensity threshold will be applied to the noiseROI volume Step 2 : n_voxel_crop = X After the threshold, an erosion of X voxels will performed on the noiseROI. -> In your case with no cropping, some voxels in your final noiseROI may still be be grey matter. I would try with at least n_voxel_crop = 1 to be sure to remove voxels with partial volume effect.

With respect to the cropping, I set it to 0 for now, see my reply to Lars above where I ask whether the eroding (or clean-up) that takes place in the segmentation step is sufficient or whether I need to do additional cropping (the second message in this issue).

About the clean-up, if you are talking about the parameter in Segement > Warping & MRF > Clean Up, I never played with it.

Best, Benoît

NinaDolfen commented 5 years ago

Thank you Benoit!

It looks better when changing the cropping from 0 to 1 (but still worse than the maps without including noise ROIS). From cropping 1 to 2 I see limited changes. (and when I set cropping to 3 than there are not enough voxels anymore). It seems that 1 seems the best choice here.

Your suggestions definitely helped, but I guess I need to tweak more settings to get it good. Did you include noise ROIS in your models before for task-related fMRI? Did it 'work' in your case?

bests,

Nina

benoitberanger commented 5 years ago

In my case, the noise ROIs approach did work in general. I have 2 different usage :

1) On pilot study, because my job is mostly setting up protocol, usually I have 2 effects :

2) On the whole cohort, I only have 1 experience in fully analyzing a study. This protocol objective was to study the conscious vs unconscious breathing brain network. The main problem is the breathing network task regressors are highly correlated to the noise induced by breathing artifacts. All the activation maps of first and second level where polluted by classic breathing artifacts. But, after using the noise ROI in WM & CSF, artifacts disappeared, mostly at the second level, and the group activation of the breathing network where clean.

My standard approach with PhysIO is to prepare 3 models with different noise regressors :

NinaDolfen commented 5 years ago

Thank you Benoit!

Can I check a few final things with you? I want to rule out that maybe there is an issue with 'how' I normalise the tissue class images.

Sorry for all the questions. I just would really like to get it to work.

Nina

NinaDolfen commented 5 years ago

Hello,

I tried a different approach for the segmentation of the tissue maps (normalising them to the dartel template (MNI space) like I normalise my structural to dartel). This did the job!

I decided to only include the CSF noise ROI (average + 2 components) because the inclusion of the WM noise ROI keeps inducing more noise in my maps.

I would like to thank again for all the input!!

Nina

benoitberanger commented 5 years ago

the tissue maps you use as input in PhysIO toolbox, are these the maps coming from the SPM segmentation step? If yes, could you specify which output exactly you use? (usually the output options are: cr.nii (for Dartel import), c.nii (native), wc.nii (normalized) or even wmc.nii (modulated))

Yes, I use or SPM12 standard Segmentation job, or the CAT12 version. I use the wc*.nii (normalized, from SPM12)

For the analyses of the full study you mention, do you remember whether you used a Dartel group template in the end to normalise your functional and structual images? If you did, how did you obtain your normalized tissue class images to input in the toolbox? I did the following: after the normalization to the DARTEL group template I segment the normalized structural and I use this output (the native images, c*.nii). So I have a segmentation step pre-normalisation (to get the necessary output for the normalisation) and an extra one after the normalisation.

No I don't use Dartel approach to build a template.

In the toolbox, do you set the forced coregistration option to yes or no? Did you ever notice whether you also get the warnings? So basicially, are your tissue images usually already coregistered with the functional or not when you go to the noise ROI extraction?

I switch off the "force coregistration" because I already performed the Realign:Coregister&Reslice wc[23].nii -> wmean.nii I don't get the warnings.

mrikasper commented 4 years ago

Dear Nina and Benoît,

what a great thread! Thank you for helping each other out, I also learnt a lot reading about the intricacies of defining the tissue probability maps.

I will close this issue soon, if there is nothing to add.

All the best, Lars