Example of a volume-based searchlight
import AnatSearchlight.searchlight as sl
import numpy as np
import nibabel as nib
import rsatoolbox as rsa
To define a searchlight we need:
a functional mask image (from the GLM) which specifies the voxels availabe for searchlight analysis.
You can define the searchlight using either:
a specific radius in mm, or
a specific number of voxels.
To specify only one of these parameters, set the other to np.inf.
# Give it an appropriate structure name: Used for generating nifti-files for output
mySearchlight = sl.SearchlightVolume('whole_brain')
# This is the mask image for the functional data space
voxel_mask = 'example_data/sub-02_ses-s1_mask.nii'
# Here we defining the searchlight parameters
mySearchlight.define(roi_img=voxel_mask,mask_img= None,maxradius=10,maxvoxels=np.inf)
mySearchlight.save('example_data/searchlight_volume.h5')
Processing center 0 of 70752
Processing center 1000 of 70752
Processing center 2000 of 70752
Processing center 3000 of 70752
Processing center 4000 of 70752
Processing center 5000 of 70752
Processing center 6000 of 70752
Processing center 7000 of 70752
Processing center 8000 of 70752
Processing center 9000 of 70752
Processing center 10000 of 70752
Processing center 11000 of 70752
Processing center 12000 of 70752
Processing center 13000 of 70752
Processing center 14000 of 70752
Processing center 15000 of 70752
Processing center 16000 of 70752
Processing center 17000 of 70752
Processing center 18000 of 70752
Processing center 19000 of 70752
Processing center 20000 of 70752
Processing center 21000 of 70752
Processing center 22000 of 70752
Processing center 23000 of 70752
Processing center 24000 of 70752
Processing center 25000 of 70752
Processing center 26000 of 70752
Processing center 27000 of 70752
Processing center 28000 of 70752
Processing center 29000 of 70752
Processing center 30000 of 70752
Processing center 31000 of 70752
Processing center 32000 of 70752
Processing center 33000 of 70752
Processing center 34000 of 70752
Processing center 35000 of 70752
Processing center 36000 of 70752
Processing center 37000 of 70752
Processing center 38000 of 70752
Processing center 39000 of 70752
Processing center 40000 of 70752
Processing center 41000 of 70752
Processing center 42000 of 70752
Processing center 43000 of 70752
Processing center 44000 of 70752
Processing center 45000 of 70752
Processing center 46000 of 70752
Processing center 47000 of 70752
Processing center 48000 of 70752
Processing center 49000 of 70752
Processing center 50000 of 70752
Processing center 51000 of 70752
Processing center 52000 of 70752
Processing center 53000 of 70752
Processing center 54000 of 70752
Processing center 55000 of 70752
Processing center 56000 of 70752
Processing center 57000 of 70752
Processing center 58000 of 70752
Processing center 59000 of 70752
Processing center 60000 of 70752
Processing center 61000 of 70752
Processing center 62000 of 70752
Processing center 63000 of 70752
Processing center 64000 of 70752
Processing center 65000 of 70752
Processing center 66000 of 70752
Processing center 67000 of 70752
Processing center 68000 of 70752
Processing center 69000 of 70752
Processing center 70000 of 70752
You can now load the pre-computed searchlight to inspect its properties.
# Here we wish to inspect the average radius and average number of voxels in each center
mySearchlight = sl.load('example_data/searchlight_volume.h5')
print(f'Average size of the searchlight: {mySearchlight.radius.mean():.1f} mm')
print(f'Average number of voxels: {mySearchlight.nvoxels.mean():.1f}')
Average size of the searchlight: 9.9 mm
Average number of voxels: 207.3
To run a searchlight we need:
an MVPA function
datafiles
The searchlight will calculate the MVPA function on data from datafiles in each center defined in mySearchlight.
# Here we use simple example that calculates the mean activity in the searchlight region.
# The MVPA functions take a `n_images x n_voxels` array as input and return either a scalar or vector as output
def mvpa_mean_function(data):
# Example of MVPA function that returns a scalar as an output argument
return np.mean(data)
def mvpa_multi_function(data):
# Example of MVPA function that returns a vector as an output argument
return np.mean(data,axis=1)
## Run the searchlight with the defined scalar MVPA function
datafiles = [f"example_data/sub-02_ses-s1_run-{r:02d}_reg-{s:02d}_beta.nii" for r in range(1,3) for s in range(2)]
results = mySearchlight.run(datafiles, mvpa_mean_function,verbose=False)
mySearchlight.data_to_nifti(results,outfilename='example_data/output_scalar.nii')
<nibabel.nifti1.Nifti1Image at 0x2846ac651f0>
## Run the searchlight with the defined scalar vector function
datafiles = [f"example_data/sub-02_ses-s1_run-{r:02d}_reg-{s:02d}_beta.nii" for r in range(1,3) for s in range(2)]
results = mySearchlight.run(datafiles, mvpa_multi_function,verbose=False)
mySearchlight.data_to_nifti(results,outfilename='example_data/output_vector.nii')
<nibabel.nifti1.Nifti1Image at 0x28464e31430>
Running representational similarity analysis (RSA) as the MVPA function:
here, we will use the rsatoolbox to run a representational similarity analysis (RSA) as the MVPA function. The first step is to prewhiten the betas using the residual mean square (resms) image from the GLM. Then we will calculate a crossnobis RDM.
# append resms to datafiles
datafiles.append('example_data/sub-02_ses-s1_resms.nii')
# define function to calculate
def mvpa_mean_rsa(data,cond_vector, partition_vector):
# 1.prewhite the betas
y = np.zeros_like(data)
for i in range(data.shape[0] - 1):
resms = data[-1,:]
resms[resms <= 0] = np.nan
y[i,:] = data[i,:] / np.sqrt(resms)
y = y[:-1,:] # remove the resms from the data
# 2.Calculate the RDM using rsatoolbox
dataset = rsa.data.Dataset(y,channel_descriptors={'channel': np.array(['vox_' + str(x) for x in range(y.shape[-1])])},obs_descriptors={'conds': cond_vector,'run': partition_vector})
rdm = rsa.rdm.calc_rdm(dataset, method='crossnobis', descriptor='conds', cv_descriptor='run')
rdm.rdm_descriptors = {'index': [0]}
# 3.get the upper triangular
rdm_vec = rdm.dissimilarities
return np.mean(rdm_vec)
To calculate the crossnobis RDM, we need to provide condition and partition vectors. The condition vector specifies the experimental condition for each beta, and the partition vector specifies the run from which each beta was estimated.
condition_vector = [1,2,1,2] # condition number for each beta
partition_vector = [1,1,2,2] # run number for each beta
result = mySearchlight.run(datafiles,mvpa_mean_rsa,function_args={'cond_vector': condition_vector, 'partition_vector': partition_vector},verbose=True)
mySearchlight.data_to_nifti(result,outfilename='example_data/output_scalar_RSA.nii')
Loading data
Volume number 1
Volume number 2
Volume number 3
Volume number 4
Volume number 5
Calculating searchlight 0 of 70752
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
Calculating searchlight 1000 of 70752
Calculating searchlight 2000 of 70752
Calculating searchlight 3000 of 70752
Calculating searchlight 4000 of 70752
Calculating searchlight 5000 of 70752
Calculating searchlight 6000 of 70752
Calculating searchlight 7000 of 70752
Calculating searchlight 8000 of 70752
Calculating searchlight 9000 of 70752
Calculating searchlight 10000 of 70752
Calculating searchlight 11000 of 70752
Calculating searchlight 12000 of 70752
Calculating searchlight 13000 of 70752
Calculating searchlight 14000 of 70752
Calculating searchlight 15000 of 70752
Calculating searchlight 16000 of 70752
Calculating searchlight 17000 of 70752
Calculating searchlight 18000 of 70752
Calculating searchlight 19000 of 70752
Calculating searchlight 20000 of 70752
Calculating searchlight 21000 of 70752
Calculating searchlight 22000 of 70752
Calculating searchlight 23000 of 70752
Calculating searchlight 24000 of 70752
Calculating searchlight 25000 of 70752
Calculating searchlight 26000 of 70752
Calculating searchlight 27000 of 70752
Calculating searchlight 28000 of 70752
Calculating searchlight 29000 of 70752
Calculating searchlight 30000 of 70752
Calculating searchlight 31000 of 70752
Calculating searchlight 32000 of 70752
Calculating searchlight 33000 of 70752
Calculating searchlight 34000 of 70752
Calculating searchlight 35000 of 70752
Calculating searchlight 36000 of 70752
Calculating searchlight 37000 of 70752
Calculating searchlight 38000 of 70752
Calculating searchlight 39000 of 70752
Calculating searchlight 40000 of 70752
Calculating searchlight 41000 of 70752
Calculating searchlight 42000 of 70752
Calculating searchlight 43000 of 70752
Calculating searchlight 44000 of 70752
Calculating searchlight 45000 of 70752
Calculating searchlight 46000 of 70752
Calculating searchlight 47000 of 70752
Calculating searchlight 48000 of 70752
Calculating searchlight 49000 of 70752
Calculating searchlight 50000 of 70752
Calculating searchlight 51000 of 70752
Calculating searchlight 52000 of 70752
Calculating searchlight 53000 of 70752
Calculating searchlight 54000 of 70752
Calculating searchlight 55000 of 70752
Calculating searchlight 56000 of 70752
Calculating searchlight 57000 of 70752
Calculating searchlight 58000 of 70752
Calculating searchlight 59000 of 70752
Calculating searchlight 60000 of 70752
Calculating searchlight 61000 of 70752
Calculating searchlight 62000 of 70752
Calculating searchlight 63000 of 70752
Calculating searchlight 64000 of 70752
Calculating searchlight 65000 of 70752
Calculating searchlight 66000 of 70752
Calculating searchlight 67000 of 70752
Calculating searchlight 68000 of 70752
Calculating searchlight 69000 of 70752
Calculating searchlight 70000 of 70752
<nibabel.nifti1.Nifti1Image at 0x2846aea58b0>