Nulls with parcellated data error

[rtang2100]

Description of issue

Hi Neuromaps Team,

Thanks for developing this wonderful tool! I have a question about spatial nulls significance testing with parcellated data. I have a self-generated array that basically includes a numerical value for each parcel (Schaefer 200), and would like to correlate it with other maps. I did not run into issue with correlations in neuromaps (my array and parcellated abagen map), but I ran into the following error when running spatial nulls:

rotated = nulls.alexander_bloch(abagen_parc, atlas='fsaverage', density='10k',
n_perm=100, seed=1234, parcellation=parc)
Traceback (most recent call last):
File "", line 1, in
File "/home/rtang/.local/lib/python3.8/site-packages/neuromaps/nulls/nulls.py", line 134, in alexander_bloch
coords, hemi = get_parcel_centroids(surfaces,
File "/home/rtang/.local/lib/python3.8/site-packages/neuromaps/nulls/spins.py", line 124, in get_parcel_centroids
for n, (parc, surf) in enumerate(zip(parcellation, surfaces)):
TypeError: 'Parcellater' object is not iterable

My 'parc' is generated with the following command:
schaefer = nntdata.fetch_schaefer2018('fsaverage')['200Parcels17Networks']
parc = Parcellater(annot_to_gifti(schaefer), 'fsaverage')

Thanks so much!
Catherine

Code of Conduct

• I agree to follow the neuromaps Code of Conduct

[rtang2100]

I actually was able to resolve this issue, however i'm running into a new error:

rotated = nulls.alexander_bloch(abagen_parc, atlas='fsaverage', density='10k',
... n_perm=100, seed=1234, parcellation=parc)
Traceback (most recent call last):
File "", line 1, in
File "/home/rtang/.local/lib/python3.8/site-packages/neuromaps/nulls/nulls.py", line 134, in alexander_bloch
coords, hemi = get_parcel_centroids(surfaces,
File "/home/rtang/.local/lib/python3.8/site-packages/neuromaps/nulls/spins.py", line 136, in get_parcel_centroids
roi = np.atleast_2d(vertices[mask].mean(axis=0))
IndexError: boolean index did not match indexed array along dimension 0; dimension is 10242 but corresponding boolean dimension is 163842

the command I used to generate abagen_pac was:

chaefer = nntdata.fetch_schaefer2018('fsaverage')['200Parcels17Networks']

parc = annot_to_gifti(schaefer)
print(parc)
(<nibabel.gifti.gifti.GiftiImage object at 0x7f6545aff1f0>, <nibabel.gifti.gifti.GiftiImage object at 0x7f6545aff490>)
newparc= parcellate.Parcellater(parc, 'fsaverage').fit()
abagen_parc = newparc.transform(abagen, 'fsaverage')
print(abagen_parc.shape)
(200,)

[justinehansen]

Hi Catherine, thanks for reaching out!
The parc input should actually be a path to the parcellation file, mapping vertices or voxels to each parcel of your atlas, defined in a specific space. Then space and den refer to the space and density of the parc file, not the annotation that you rotate. This isn't so clear in the documentation so that's our bad, sorry about that!

Here's what you should do instead:

# fetch the parcellation file; in this case this is in fsaverage 164k space
schaefer = nntdata.fetch_schaefer2018('fsaverage')['200Parcels17Networks']

# make gifti
schaefer = annot_to_gifti(schaefer)

# get nulls
rotated = nulls.alexander_bloch(data=abagen_parc, atlas='fsaverage', density='164k', n_perm=1000, seed=1234, parcellation=schaefer)

Hope that helps! Let me know if you have any other questions 🙂

Best,
Justine

[rtang2100]

Hi Justine,

Thanks a lot for prompt response! The issue is resolved now. ;)
One last question: is there a way for me to get the distribution of the spatial nulls? I tried out print (rotated), but I'm not positive that it gave me the correct info as it included values that are >abs(1). I'd like to plot it out like Fig. 4 in the neuromaps paper.

Thanks,
Catherine

[justinehansen]

Hi Catherine,

Unfortunately neuromaps doesn't return the nulls but I completely agree that we should 😅 I also like to plot the null distributions and so far I always do it manually. For now, while it isn't yet implemented in neuromaps, here's how I do it in my own work:

# get data arrays
mymap_arr = load_data(mymap)
targetmap_arr = load_data(targetmap)
rho = pearsonr(mymap_arr, targetmap_arr)[0]  # same output as from `compare_images(mymap, targetmap, metric='pearsonr')

# get rotated data
nspins = 10000  # number of rotations you want to do
mymap_rot = nulls.alexander_bloch(data=mymap, atlas='fsaverage', density='164k', n_perm=nspins, seed=1234, parcellation=schaefer)

# run spin test "manually" to save null distributions for plotting
null_rho = np.zeros((nspins, ))
for i in range(nspins):
    null_rho[i] = pearsonr(mymap_rot[:, i], targetmap_arr)[0]
    
# calculate p-value
pspin = (1 + sum(abs(null_rho - np.mean(null_rho)) > abs(rho - np.mean(null_rho)))) / (nspins + 1)

I'll be implementing this in neuromaps soon so we can fetch the distribution straight from the compare_images function! By the way, the rotated output is your parcellated data but randomly shuffled in a way that preserves spatial autocorrelation (ie the data is being spun around and the output is "rotated"). However, if you set data=None in the function (rotated = nulls.alexander_bloch(data=None, atlas='fsaverage', density='164k', n_perm=nspins, seed=1234, parcellation=schaefer)), it will return the indices to rotate your data, instead of the rotated data itself. This is convenient because these indices can be reused to rotate other maps that are defined in the same Schaefer-200 parcellation.

Let me know if you have any other questions!

Justine

[rtang2100]

Hi Justine,

Thanks so much for these helpful code and information! Really appreciated! :)

Best,
Catherine

[justinehansen]

Hi @rtang2100 just letting you know that I've just merged a PR that makes it possible to return the null distribution from compare_images(). Hopefully that makes things easier!

Best,
Justine