weixin_39809540
weixin_39809540
2021-01-01 11:45

FIX: Revise the reproducibility of CompCor masks

This PR revises the implementation of CompCor masks in an attempt to make it closer to the original proposal and at the same time address some recurring problems of fMRIPrep's tCompCor implementation.

Finally, with the more careful resampling of prior knowledge from the anatomical scan, this refactor should also make the aCompCor components more run-to-run repeatable.

aCompCor

The massaging of CompCor masks is now done in anatomical space where it is more precise, and a careful resampling to BOLD space follows.

The implementation deviates from Behzadi et al. Their original implementation thresholded the CSF and the WM partial-volume masks at 0.99 (i.e., 99% of the voxel volume is filled with a particular tissue), and then binary eroded that 2 voxels:

Anatomical data were segmented into gray matter, white matter, and CSF partial volume maps using the FAST algorithm available in the FSL software package (Smith et al., 2004). Tissue partial volume maps were linearly interpolated to the resolution of the functional data series using AFNI (Cox, 1996). In order to form white matter ROIs, the white matter partial volume maps were thresholded at a partial volume fraction of 0.99 and then eroded by two voxels in each direction to further minimize partial voluming with gray matter. CSF voxels were determined by first thresholding the CSF partial volume maps at 0.99 and then applying a threedimensional nearest neighbor criteria to minimize multiple tissue partial voluming. Since CSF regions are typically small compared to white matter regions mask, erosion was not applied.

This particular procedure is not generalizable to BOLD data with different voxel zooms as the mathematical morphology operations will be scaled by those. Also, from reading the excerpt above and the tCompCor description, I () believe that they always operated slice-wise given the large slice-thickness of their functional data.

Instead, fMRIPrep's implementation deviates from Behzadi's implementation on two aspects:

  • the masks are prepared in high-resolution, anatomical space and then projected into BOLD space; and,
  • instead of using binary erosion, a dilated GM map is generated -- thresholding the corresponding PV map at 0.05 (i.e., pixels containing at least 5% of GM tissue) and then subtracting that map from the CSF, WM and CSF+WM (combined) masks. This should be equivalent to eroding the masks, except that the erosion only happens at direct interfaces with GM.

When the probseg maps provene from FreeSurfer's recon-all (i.e., they are discrete), binary maps are transformed into some sort of partial volume maps by means of a Gaussian smoothing filter with sigma adjusted by the size of the BOLD data.

tCompCor

In the case of tCompCor, this commit removes the heavy erosion of the brain mask because 1) that wasn't part of the original proposal by Behzadi et al., and 2) the erosion was the potential source of errors from numpy complaining that it can't take from an empty axis of an array.

Based on these results, we chose a 2% threshold (∼20–30 voxels per slice) as a reasonable empirical threshold that effectively identified voxels with the highest fractional variance of physiological noise.

Although they do the calculation slice-wise, this commit rolls tCompCor back to calculate the 2% threshold on the whole-brain mask.

I would appreciate a deep review from & to ensure my interpretation of Behzadi's paper is correct. I'm also tagging for a more top-level check of the changes. Further feedback from others is of course very much appreciated.

Resolves: #2129. References: #2052.

该提问来源于开源项目:nipreps/fmriprep

  • 点赞
  • 写回答
  • 关注问题
  • 收藏
  • 复制链接分享
  • 邀请回答

19条回答

  • weixin_39663729 weixin_39663729 4月前

    Pushed your previous commit to upstream/compcormasks to generate reports. The good news is: my changes didn't cause the bad-looking reports.

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    The masks in reports look pretty bad.

    You mean for ds054 and ds210 (i.e., --fs-no-reconall), right?

    I have mostly tested with recon-all - that's why I haven't lent much attention to ds054 and ds210. That said, I'm surprised that the combined masks look so bad.

    To recapitulate:

    Before this PR we were probably representing the combined mask (as stated in the caption of the reportlet), and it never ever contained CSF. I think our CSF mask was so exclusive that very few voxels were in it (and hence, the empty take errors from numpy). - I say this just in case you think ds005's reports look awful too.

    Now, the masks with recon-all look pretty good. I've checked all of them after running ds114, and has also generated a good bunch (which looked good to me and will soon assess after his quals).

    Then, we have the --fs-no-reconall path. I know I have run ds054 without reconall at full resolution and that gave me sound results. This might be a result of the very big pixels of the WM mask with the heavily subsampled dataset. Let me run a bunch of dataset. Perhaps, I should start with ds114 (and --fs-no-reconall), and see how different the masks are w.r.t. the default.

    点赞 评论 复制链接分享
  • weixin_39663729 weixin_39663729 4月前

    You're right, I didn't look at ds005 because it takes so much longer to complete.

    Do the voxels near the edge of the brain mask occur in higher resolution images? Are they problematic? I assume they're classified as CSF, but I would worry you get GM signal bleed-over (possibly resulting in the high variance for tCompCor).

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    Do the voxels near the edge of the brain mask occur in higher resolution images? Are they problematic? I assume they're classified as CSF, but I would worry you get GM signal bleed-over (possibly resulting in the high variance for tCompCor).

    Are you referring to pink (aCompCor) or blue (tCompCor)? I guess they are not pink because the new implementation is quite careful - if they are selected, the certainty of them being CSF is very high.

    Regarding tCompCor, this is expected (which does not mean good). I particularly asked about this, and he didn't seem too concerned. Please note that the original Behzadi paper DID NOT erode the brain mask (as we were doing), and also, their top percentile was calculated slice-wise (with veeery thick slices). My impression is that the original implementation of tCompCor wasn't great because of this.

    We may want to come back to the previous implementation, but I definitely see the new tCompCor masks picking up voxels in the brain stem and external parts of the commissures - which is a good thing, I believe.

    点赞 评论 复制链接分享
  • weixin_39663729 weixin_39663729 4月前

    Are you referring to pink (aCompCor) or blue (tCompCor)?

    They show up in both.

    Fair enough.

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    Okay, running ds054 and ds210 at their original resolution with --fs-no-reconall. Will report back when done.

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    Okay, running ds054 and ds210 at their original resolution with --fs-no-reconall. Will report back when done.

    Masks look okay (better with FS, but that's the drill). I'm building a new singularity image (which required cutting an RC of sMRIPrep) to test these masks including the fix from Chris to make the reportlet show the combined mask instead of the wm mask.

    As soon as I have both datasets run with the latest image I'll ask for his trained eye.

    Meanwhile, please have a look into #2245 - I think it would be worth pushing into 20.1.x to be able to do the correlation w/global-signal test.

    点赞 评论 复制链接分享
  • weixin_39663729 weixin_39663729 4月前

    What's the status of this? I think we're very close to a release, so let me know if I can help on anything.

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    My plan is to quickly make a PR to add the CompCor time series to the bar plot of correlations with the GS (and remove the components from there). That if thinks it's a good idea.

    Worst case scenario, the plots will become more readable.

    On Wed, Sep 2, 2020, 14:57 Chris Markiewicz wrote:

    https://github.com/oesteban https://github.com/rciric What's the status of this? I think we're very close to a release, so let me know if I can help on anything.

    — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/poldracklab/fmriprep/pull/2130#issuecomment-685717634, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAESDRSHLSWESN4TUY4AXRDSDY6L7ANCNFSM4NGMUUGQ .

    点赞 评论 复制链接分享
  • weixin_39663729 weixin_39663729 4月前

    Is this PR ready to merge, then?

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    Yes, I think that's safe

    On Wed, Sep 2, 2020, 17:06 Chris Markiewicz wrote:

    Is this PR ready to merge, then?

    — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/poldracklab/fmriprep/pull/2130#issuecomment-685797903, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAESDRXJYUT22W32H2R5QNTSDZNOPANCNFSM4NGMUUGQ .

    点赞 评论 复制链接分享
  • weixin_39888180 weixin_39888180 4月前

    Hello,

    Apologies for the late response. I agree that morphological erosion using an isotropic fixed-voxel-size structuring element would generalize poorly to new resolutions. Reading the Behzadi paper again, I’m somewhat confused about the decision to warp the data into functional space before applying the erosion, and the motivation for doing so is unclear to me. I completely agree with the decision to perform the erosion in anatomical space, and the adaptive kernel size is a great idea.

    Eroding using a low-threshold GM mask intuitively sounds like an improvement over conventional methods, but I do think that we should ideally examine the correlations of WM and CSF time series with whole-brain global signal when using the new approach to determine whether there are any major changes from the conventional approach. In particular, we would expect that limiting partial volume effects keeps the tissue-specific correlation with global signal reasonably low. This is of particular concern for those researchers who believe strongly that global mean signal should not be discarded (although I do not find myself in that camp).

    My recommendation for tCompCor is similar: the proposed approach sounds very reasonable, although I think it wouldn't hurt to more clearly characterise the noise ROI and the component time series correlations at some point. tCompCor appears to be used with somewhat less frequency than aCompCor.

    To verify, are the new aCompCor masks also going to define the regions over which the (global) mean WM and CSF signals are extracted?

    The literature is not entirely uniform in use of anatomical CompCor approaches. In addition to the original method from Behzadi and colleagues, a similar aCompCor methodology proposed by Muschelli and colleagues is used with considerable frequency. In particular, the aCompCor evaluated by our group (and by Parkes et al. to my knowledge) more closely matched the Muschelli approach. By computing separate decompositions for each compartment mask and for the union mask, fMRIPrep already supports the most significant difference between the approaches. It's not, however, immediately clear to me whether Muschelli and colleagues performed their mask erosions in native anatomical space or in a standard space, since the maps they used were provided by a combined segmentation-normalisation procedure from SPM:

    Subject-specific tissue probability maps generated during SPM’s unified segmentation-normalization process were restricted using a 99% probability threshold. To further reduce the risk of capturing signals of interest from adjacent grey matter voxels, the resulting WM mask was eroded using MATLAB’s imerode function, and the CSF mask was constrained to areas within the ALVIN mask of the ventricles (Kempton et al., 2011).

    As a potential aid for evaluating the masks we obtain, I'm including here for reference the subject-wise voxel inclusion frequency map for WM and CSF masks in the Muschelli study. image

    (from Muschelli et al., 2014)

    Also potentially of note, Power and colleagues performed their mask erosion using FreeSurfer. I'm not sure whether this erosion was performed in volumetric or surface space, and I only have a passing familiarity with FreeSurfer. In any case, there is substantial precedent for eroding masks in subject anatomical space. This approach was most similar to what our group did in the benchmarking paper, although we performed our erosions morphologically in volumetric space.

    Each subject’s MP-RAGE underwent segmentation by FreeSurfer (version 5.0). Segments corresponding to the lateral ventricles, the white matter, the gray matter, and all within-brain voxels were extracted. 5 versions of the white matter and ventricular masks were formed: uneroded, and erosions 1–4 times. All masks were resliced to 3 mm isotropic dimensions to match the BOLD data. Maximal erosion (4×) masks were preferred; if maximal erosion yielded empty masks, lesser erosions were progressively considered until a mask with qualifying voxels was found. This relaxation occurred infrequently for white matter masks. Erosions of 1 were often used for CSF masks. These masks were used to extract nuisance signals, to calculate DV and SD within the masks, and to define the voxel timeseries displayed in Part I.

    The description in a later work suggests to me that the erosion was likely also in volumetric space, although I again cannot say with confidence. This study also includes an example of eroded masks (corresponding to the deepest WM colour below, although I seem to recall that many of the example masks I saw from the WashU group were restricted to deeper voxels still.)

    image

    (from Power et al., 2017)

    Overall, I agree with the premises of this revision, although I believe that further practical characterisation of fMRIPrep's CompCor procedure, potentially in the form of empirical correlations and voxel inclusion maps, would be welcome at some future stage.

    点赞 评论 复制链接分享
  • weixin_39663729 weixin_39663729 4月前

    I don't think the FreeSurfer method described could apply to the surface, as there is no representation of white matter or CSF on the surface. It sounds like they used the aseg segmentation.

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    Hi , happy to do this:

    we should ideally examine the correlations of WM and CSF time series with whole-brain global signal when using the new approach to determine whether there are any major changes from the conventional approach. In particular, we would expect that limiting partial volume effects keeps the tissue-specific correlation with global signal reasonably low.

    If I understand correctly and to say that this PR improves the situation, we should calculate the correlation of global signal with the time-series extracted from the different CompCor masks. Then, if these correlations are not significantly larger than those before this PR, we can say these new masks at least don't make worse masks for compcor. Is that right?

    点赞 评论 复制链接分享
  • weixin_39663729 weixin_39663729 4月前

    The masks in reports look pretty bad.

    点赞 评论 复制链接分享
  • weixin_39748773 weixin_39748773 4月前

    Thank your for raising your pull request.

    Some of the fMRIPRep maintainers will review your changes as soon as time permits. I'm attaching below a Review Checklist for the reviewers, to check off during the review.

    PR Review

    Please check off boxes as applicable, and elaborate in comments below. Your review is not limited to these topics, as described in the reviewer guide

    • [ ] As the reviewer I confirm that there are no conflicts of interest for me to review this work.

    Please check what applies in the following aspects of the PR:

    Code documentation

    • [ ] New functions and modules are documented following the guidelines.
    • [ ] Existing code required amendments in documentation and has been updated.
    • [ ] Sufficient inline comments ensure readability by other contributors in the long term.

    Documentation site

    • [ ] The modifications to fMRIPrep in this PR do not require extra documentation. This is a common situation when a bug fix that does not change the code base substantially is submitted.
    • [ ] This PR requires documentation. The corresponding documentation will be submitted as an additional PR in the future.
    • [ ] This PR requires extra documentation, and will be included within this PR. Please check the boxes that apply best:
    • [ ] An existing feature is substantially modified, changing the interface and/or logic of a module.
    • [ ] A new feature is being added, therefore full documentation of it will be included
    • [ ] This PR addresses an error in existing documentation.
    • [ ] Yes, all expected sections of documentation were generated by the CI build, and look as expected

    Tests

    • [ ] The new functionalities in this PR are covered by the existing tests
    • [ ] New test batteries have been included with this PR

    Data

    • [ ] This PR does not require any additional data to be uploaded to OSF.
    • [ ] This PR requires additional data for tests.

    Dependencies: smriprep

    • [ ] This PR does not require any update on smriprep; or
    • [ ] smriprep has correctly been pinned.

    Dependencies: niworkflows

    • [ ] This PR does not require any update on niworkflows; or
    • [ ] niworkflows has correctly been pinned.

    Dependencies: sdcflows

    • [ ] This PR does not require any update on sdcflows; or
    • [ ] sdcflows has correctly been pinned.

    Dependencies: Nipype

    • [ ] This PR does not require any update on nipype; or
    • [ ] nipype has correctly been pinned.

    Dependencies: other

    • [ ] This PR does not require any update on other dependencies; or
    • [ ] other dependencies have correctly been pinned.

    Reports generated within CI tests

    • [ ] Yes, I have checked the reports of ds005 and they are not modified or the changes are as expected
    • [ ] Yes, I have checked the reports of ds054 and they are not modified or the changes are as expected
    • [ ] Yes, I have checked the reports of ds010 and they are not modified or the changes are as expected
    点赞 评论 复制链接分享
  • weixin_39821874 weixin_39821874 4月前

    Hello ! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

    There are currently no PEP 8 issues detected in this Pull Request. Cheers! :beers:

    Comment last updated at 2020-08-28 08:00:36 UTC
    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    Ready for you

    点赞 评论 复制链接分享
  • weixin_39809540 weixin_39809540 4月前

    The horrible contours on our subsampled dataset will be addressed by nipreps/niworkflows#544.

    点赞 评论 复制链接分享

相关推荐