I am trying to recreate analyses included in this paper (Grotheer et al., 2019) but with Python/CLI tools. So my goal is to do the following:
Find functional ROIs that correspond to task activation for particular tasks (done)
a) Used fMRIPrep and SPM
Construct whole brain tractography seeded from the gray matter white matter interface (done)
b) QSIPrep for preprocessing and mrtrix3 for tractography
Segment the tractography into candidate fiber bundles (done)
c) DIPY Recobundles
Find the GMWMI adjacent to the function ROIs (???)
Find streamlines within the candidate fiber bundles that pass through the GMWMI for a given functional ROI (???)
Use these streamlines for analyses.
I am stuck on steps 4 and 5. For step 4, I have my GMWMI map used for tractography and the functional ROIs, but I don’t know how to find the portion of GMWMI that corresponds to a given functional ROI. For step 5, I have all of my candidate fiber tracts derived from the whole brain tractogram, but do not know how to find only the streamlines that pass through a certain ROI (which would be the output of step 4).
Any guidance would be greatly appreciated, and I’d be happy to provide more details on my implementation as needed.
For step 4, when you say “Find the GMWMI adjacent to the function ROIs” does the term “adjacent” mean nearby or does it mean overlapping?
In the case that you mean overlapping, then this should be a fairly straightforward “AND” operation between the two binarized masks corresponding to the GMWMI and the functional ROI. I think this nilearn function should do the trick : nilearn.masking.intersect_masks
I think you will find that your issue with #5 is a little underspecified. If you’re only using one ROI as a segmentation criteria (GM based or not) you’re typically going to end up with A LOT of streamlines going all over the place (the brain is a fairly densely connected thing, after all). As such, you’ll likely need an additional set of criteria to narrow down your candidate streamlines. The way you wish to narrow this down may depend on your overarching research question.
Insofar as just using a ROI as a segmentation criteria, I’m not sure if there’s any standard code for this in python. There is some in matlab though with the vistasoft package. I’ve hacked together a python equivalent of this here. Also there may be something comparable behind the scenes in pyAFQ but I’m not sure.
My response here is a bit disorganized, but let me know if there’s anything else I can do to help.
In regards to step 4: I suppose that will depend on the location of the functional ROIs. If there is an overlap between the GM and GMWMI then that would be great, but since that isn’t necessarily the case, it would have to be the nearest portion of the GMWMI.
Step 5: definitely agree with what you said. My approach was to only consider fibers from a few bundles (AF, SLF, IFOF, ILF) and find fibers within these bundles that intersect each ROI. I’ll look into the links you posted; should be very helpful.
A potential solution to the issue with #4 is to use scipy’s dilate function. If you’re trying to find the nearest GM area, then I think inflating/dilating the fROI should work (and then using the intersect function).
Sounds like a reasonable solution, thanks! Looks like MrTrix has a function that can do the track intersect tckedit full_bundle.tck full_bundle_intersect.tck -include ROI_to_Intersect.mif. However, it is unclear whether this ROI_to_Intersect has to be GMWMI or original GM, since some functions in MrTrix (like tck2connectome) already do a radial search at the end of streamlines to decide which ROI to assign to, which may serve the same purpose as dilating a fROI.