[BUG] Index error when parcellating using Schaefer1000Parcels7Networks

[ashlea-segal]

Issue summary

Thank you for this excellent toolbox!

I am trying to parcellate data using the Schaefer1000Parcels7Network parcellation and get the following error message:

It seems to work for coarser resolutions so I suspect the error occurs because the Schaefer 1000 Parcel parcellation is missing two ROIs?

Thank you so much!

Detailed issue description

No response

Steps to reproduce issue

schaefer_fslr = datasets.fetch_schaefer2018(version='fslr32k')
giis = neuromaps.images.dlabel_to_gifti(schaefer_fslr['1000Parcels7Networks'])
parc = Parcellater(giis, 'fslr')
cmro2_fslr = neuromaps.datasets.fetch_annotation(desc='cmr02')
cmro2_fslr = parc.fit_transform(cmro2_fslr, 'fslr')

Software version

No response

What operating system were you using when you encountered this issue?

• macOS
• Windows
• Linux

Code of Conduct

• I agree to follow the neuromaps Code of Conduct

[VinceBaz]

Hi Ashlea! You are correct: you get this error message because some parcels are missing from the Schaefer 1000 parcellation :)

The error message was unfortunately quite vague, so I'll try to make it clearer!

Also, as they explained, the parcellation misses a parcel as a result of the fsaverage-to-fslr transformation used to generate the fslr version of the parcellation. Issues like this (i.e. some parcels disappearing) are bound to happen when high-resolution parcellations are transformed from one coordinate system to another. I would suggest avoiding transforming parcellation atlases as much as possible and transforming the brain maps instead. This way, we don't lose any parcel :)

For you specific example, for instance, I would use the fsaverage6 version of the Schaefer 1000 parcellation atlas instead. The fit_transform function would then transform the CMRO2 map to the fsaverage coordinate system before parcellating it. It would basically look as follows:

schaefer_fs6 = datasets.fetch_schaefer2018(version='fsaverage6')
giis = neuromaps.images.annot_to_gifti(schaefer_fs6['1000Parcels7Networks'])
giis = neuromaps.images.relabel_gifti(giis)
parc = Parcellater(giis, 'fsaverage', resampling_target='parcellation')
cmro2_fslr = neuromaps.datasets.fetch_annotation(desc='cmr02')
cmro2_parc = parc.fit_transform(cmro2_fslr, 'fslr')

(relabel_gifti is not required in my example, but I personally like to call it anytime I parcellate a brain map, to make sure that the labels are consecutive!)

This should work! :)

Vincent

[ashlea-segal]

Excellent. Thank you so much, Vincent!

I plan to use multiple brain maps from neuromaps across a variety of coordinate space. Where possible do you recommend parcellating the brain map from its original coordinate space, or transforming all brain maps to a common coordinate space first?

As a side - I ran both versions of the code above using coarser parcellations (e.g. Schaefer300Parcels7Networks) - perhaps naively, I expected there to be a linear relationship between the parcellated data generated using fsaverage6 Schaefer atlas and the parcellated data using fslr Schaefer atlas. For the most part, this looks to be the case, but not for all ROIs (see example below). Have you encountered this problem before? I checked that the Schaefer atlases in both spaces are similar, hence I do think that is causing the problem.

Again, really appreciate this wonderful toolbox and your very prompt response :)

[ashlea-segal]

oops sorry - I didn't mean to close the issue!

[VinceBaz]

I am glad you like the toolbox 😄. Here are my answers to your questions:

original vs. common coordinate space question:

I think that in theory both options are fine, especially since you'll parcellate the data. Personally however, I think I would prefer doing as few transformation as possible. So I would probably use the original coordinate spaces!

fs6 vs. fslr parcellation problem

I ran into the exact same problem yesterday while testing the code I gave you! And I found the issue:

Basically, the fit_transform method works by calling the vertices_to_parcels function, which computes the average of the brain map's vertices within each parcel, but excludes vertices that are NaN. Unfortunately, it does not ignore vertices that are 0: (e.g. the medial wall in the CMRO2 brain map)

Unfortunately, the fsaverage6 Schaefer parcellation does not perfectly match the boundary of the medial wall: you can look at the first figure below (here I use the Schaefer 800 version) and see that some parcels overlap with the medial wall (see the region I circled in black).

As a result, when the vertices_to_parcels function computes the means within each parcel, it considers some vertices that are within the medial wall, which have values of 0. So the mean of these parcels bordering the medial wall for the "fsaverage6" is much smaller than the mean of the same parcels in the "fslr" version, which do not overlap with the medial wall. You can see the difference between the two parcellated brain maps below ("resample to data" = fslr version and "resampled to parcellation" = fsaverage6 version)!

This is quite annoying and I am so sorry about this issue! Basically the "fslr" version is the most accurate one. For now, to fix this issue for the fsaverage6 version, what you'd have to do is manually set the medial wall to NaN values in the CMRO2 brain map (since this is what the vertices_to_parcels function ignores by default! It would look like this:

from neuromaps.resampling import resample_images
from neuromaps.parcellate import _gifti_to_array
from neuromaps.nulls.spins import vertices_to_parcels
import numpy as np

# Fetch Schaefer parcellation
schaefer_fs6 = fetch_schaefer2018(version='fsaverage6')['1000Parcels7Networks']
schaefer_fs6 = annot_to_gifti(schaefer_fs6)

# Fetch CMRO2 map
cmro2_fslr_164k = fetch_annotation(desc='cmr02')

# Resample CMRO2 to fsaverage space
cmro2_fs6, _ = resample_images(
    cmro2_fslr_164k, schaefer_fs6, 'fslr', 'fsaverage', resampling='transform_to_trg')

# Manually set values of `0` to `NaN`
cmro2_array = _gifti_to_array(cmro2_fs6)
cmro2_array[cmro2_array == 0] = np.nan

# Manually call `vertices_to_parcels` function
cmro2_parc_fs6_nan = vertices_to_parcels(cmro2_array, schaefer_fs6)

It's quite unfortunate that this has to be done manually. I'll try to make a Pull Request in the next couple of days that should fix this issue 😄. I'm thinking of adding a parameter to the vertices_to_parcels function that asks the user which values should be considered as the background in the brain map that we want to parcellate!

Hopefully my explanation makes sense, and sorry again for this issue!

Vincent

[VinceBaz]

Hey Ashlea, I opened a Pull Request that should hopefully fix the issue that I described above 😄 . Once I merge this PR, we will now be able to use the ignore_background_data parameter in the transform method of the Parcellater class to ignore vertices in the medial wall that have values of 0.

Now, by running the following lines of codes, we should get a result very similar to the one we obtain when resampling to "data" (i.e. with the fslr version of the parcellation):

schaefer_fs6 = datasets.fetch_schaefer2018(version='fsaverage6')
giis = neuromaps.images.annot_to_gifti(schaefer_fs6['800Parcels7Networks'])
giis = neuromaps.images.relabel_gifti(giis)
parc = Parcellater(giis, 'fsaverage', resampling_target='parcellation')
cmro2_fslr = neuromaps.datasets.fetch_annotation(desc='cmr02')
cmro2_parc = parc.fit_transform(cmro2_fslr, 'fslr', ignore_background_data=True)

Hopefully everything should now be fixed!! Thank you so much for sharing this bug and feel free to let us know if you have other questions or issues 😄

Vincent

[ashlea-segal]

This is wonderful. Thank you so much for such a detailed, very thorough, and incredibly speedy response, Vincent. Very much appreciated! :)

[VinceBaz]

Hey @ashlea, I merged the PR and closed this issue. You can always let us know if you have any other question 😄