[BUG] Abnormal vertex reassignment in neuromaps.nulls.alexander_bloch

[alyssadai]

Issue summary

(Posting here for traceability)
Vertex values of large swathes of cortex appear to not be reassigned in null maps generated by the neuromaps.nulls.alexander_bloch method, resulting in abnormal proportions of repeated values.

Detailed issue description

Hi neuromaps team, thank you for this wonderful tool!

While attempting to use neuromaps.nulls.alexander_bloch() to spin a map of vertex-parcel assignments for a data-driven cortical parcellation (e.g. where the data is a 81924x1 array of integer parcel labels from 1-8), I noticed an issue with the generated null maps upon plotting them on my CIVET cortical surface. Certain vertices appear to be reassigned significantly more than others, resulting in multiple parcel labels being lost in the 'spun' maps.

I was unable to reproduce this behaviour using the original Alexander-Bloch matlab spin-test toolbox (the null maps from the o.g. implementation look normal), so this appears to be a neuromaps-specific issue (or an issue of how I am using it), or possibly a bug introduced by one of the methods called behind the scenes (netneurotools.stats.gen_spinsamples()?)...

The specific data I was trying to spin can be found here (in the code below, lh_nmf-k8_allcomps.txt = lh_nmf-resid-zvert_k8_clustering.txt, and rh_nmf-k8_allcomps.txt = rh_nmf-resid-zvert_k8_clustering.txt), and the surfaces I used can be found here.

Steps to reproduce issue

import neuromaps
import numpy as np     # v1.19.2

data_left = np.loadtxt('lh_nmf-resid-zvert_k8_clustering.txt', dtype='int32')
data_right = np.loadtxt('rh_nmf-resid-zvert_k8_clustering.txt', dtype='int32')
data_arr = np.concatenate((parc_left,parc_right))
surface = ('surfaces/PS_N266_lh_average_midsurface.surf.gii', 
           'surfaces/PS_N266_rh_average_midsurface.surf.gii')

nperm=100
rotated_arr = neuromaps.nulls.alexander_bloch(data_arr, surfaces=surface, n_perm=nperm, seed=528)

# checking for presence of the 8 parcel label values from the original data in the first 3 permutations
print(np.unique(rotated_arr[:,0])) # out: [0 1 2 3 4 6 8]
print(np.unique(rotated_arr[:,1])) # out: [0 1 2 4 5 6 7 8]
print(np.unique(rotated_arr[:,2])) # out: [0 1 3 4 6 8]

# plotting some example rotations
def return_gii_tuple(surfLR): # helper function to create giftis of null maps
    surfL = surfLR[:int(len(surfLR)/2)]
    surfR = surfLR[int(len(surfLR)/2):]
    img = (neuromaps.images.construct_shape_gii(surfL,intent = 'NIFTI_INTENT_LABEL'),
           neuromaps.images.construct_shape_gii(surfR,intent = 'NIFTI_INTENT_LABEL'))
    return(img)

for i in range(3):
    plot = neuromaps.plotting.plot_surf_template(return_gii_tuple(rotated_arr[:,i]), template='civet', 
                                                 density='41k', surf='inflated', cmap='Spectral', 
                                                 mask_medial=False, avg_method='median')

Software version

0.0.3

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

[alyssadai]

Update: I didn't think to test this before, but I just noticed that if I use the built-in parameters for specifying atlas and density in neuromaps.nulls.alexander_bloch() instead of supplying my own CIVET-derived surfaces (in .gii format) using the surfaces parameter as done above, I no longer seem to encounter the issue. I would have thought that both surface specifications should have the same density (just a different geometry) so I wonder if this indicates a problem arising at the stage of conversion between CIVET .obj surfaces and .gii (also done using neuromaps functions, see second code snippet below), or how neuromaps is reading the bespoke surfaces. I've uploaded my surfaces in both formats here.

rotated_arr = neuromaps.nulls.alexander_bloch(parcellation_arr, atlas='civet', density='41k', n_perm=nperm, seed=528)

# commands used to convert my .obj surfaces to .gii for use in neuromaps
neuromaps.images.obj_to_gifti('surfaces/PS_N266_lh_average_midsurface.obj')
neuromaps.images.obj_to_gifti('surfaces/PS_N266_rh_average_midsurface.obj')