{ "cells": [ { "cell_type": "markdown", "id": "5ce97562", "metadata": {}, "source": [ "## Example of a volume-based searchlight" ] }, { "cell_type": "code", "execution_count": 1, "id": "d9b85a3b", "metadata": {}, "outputs": [], "source": [ "import AnatSearchlight.searchlight as sl\n", "import numpy as np\n", "import nibabel as nib\n", "import rsatoolbox as rsa" ] }, { "cell_type": "markdown", "id": "620c36fa", "metadata": {}, "source": [ "To define a searchlight we need:\n", "\n", "* a functional mask image (from the GLM) which specifies the voxels availabe for searchlight analysis.\n", "\n", "You can define the searchlight using either:\n", "* a specific radius in mm, or \n", "* a specific number of voxels. \n", "\n", "To specify only one of these parameters, set the other to `np.inf`. \n" ] }, { "cell_type": "code", "execution_count": 2, "id": "bcc6ef85", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Processing center 0 of 70752\n", "Processing center 1000 of 70752\n", "Processing center 2000 of 70752\n", "Processing center 3000 of 70752\n", "Processing center 4000 of 70752\n", "Processing center 5000 of 70752\n", "Processing center 6000 of 70752\n", "Processing center 7000 of 70752\n", "Processing center 8000 of 70752\n", "Processing center 9000 of 70752\n", "Processing center 10000 of 70752\n", "Processing center 11000 of 70752\n", "Processing center 12000 of 70752\n", "Processing center 13000 of 70752\n", "Processing center 14000 of 70752\n", "Processing center 15000 of 70752\n", "Processing center 16000 of 70752\n", "Processing center 17000 of 70752\n", "Processing center 18000 of 70752\n", "Processing center 19000 of 70752\n", "Processing center 20000 of 70752\n", "Processing center 21000 of 70752\n", "Processing center 22000 of 70752\n", "Processing center 23000 of 70752\n", "Processing center 24000 of 70752\n", "Processing center 25000 of 70752\n", "Processing center 26000 of 70752\n", "Processing center 27000 of 70752\n", "Processing center 28000 of 70752\n", "Processing center 29000 of 70752\n", "Processing center 30000 of 70752\n", "Processing center 31000 of 70752\n", "Processing center 32000 of 70752\n", "Processing center 33000 of 70752\n", "Processing center 34000 of 70752\n", "Processing center 35000 of 70752\n", "Processing center 36000 of 70752\n", "Processing center 37000 of 70752\n", "Processing center 38000 of 70752\n", "Processing center 39000 of 70752\n", "Processing center 40000 of 70752\n", "Processing center 41000 of 70752\n", "Processing center 42000 of 70752\n", "Processing center 43000 of 70752\n", "Processing center 44000 of 70752\n", "Processing center 45000 of 70752\n", "Processing center 46000 of 70752\n", "Processing center 47000 of 70752\n", "Processing center 48000 of 70752\n", "Processing center 49000 of 70752\n", "Processing center 50000 of 70752\n", "Processing center 51000 of 70752\n", "Processing center 52000 of 70752\n", "Processing center 53000 of 70752\n", "Processing center 54000 of 70752\n", "Processing center 55000 of 70752\n", "Processing center 56000 of 70752\n", "Processing center 57000 of 70752\n", "Processing center 58000 of 70752\n", "Processing center 59000 of 70752\n", "Processing center 60000 of 70752\n", "Processing center 61000 of 70752\n", "Processing center 62000 of 70752\n", "Processing center 63000 of 70752\n", "Processing center 64000 of 70752\n", "Processing center 65000 of 70752\n", "Processing center 66000 of 70752\n", "Processing center 67000 of 70752\n", "Processing center 68000 of 70752\n", "Processing center 69000 of 70752\n", "Processing center 70000 of 70752\n" ] } ], "source": [ "# Give it an appropriate structure name: Used for generating nifti-files for output\n", "mySearchlight = sl.SearchlightVolume('whole_brain')\n", "# This is the mask image for the functional data space\n", "voxel_mask = 'example_data/sub-02_ses-s1_mask.nii'\n", "\n", "# Here we defining the searchlight parameters\n", "mySearchlight.define(roi_img=voxel_mask,mask_img= None,maxradius=10,maxvoxels=np.inf)\n", "mySearchlight.save('example_data/searchlight_volume.h5')" ] }, { "cell_type": "markdown", "id": "fa90c25c", "metadata": {}, "source": [ "You can now load the pre-computed searchlight to inspect its properties." ] }, { "cell_type": "code", "execution_count": 3, "id": "c8074c67", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Average size of the searchlight: 9.9 mm\n", "Average number of voxels: 207.3\n" ] } ], "source": [ "# Here we wish to inspect the average radius and average number of voxels in each center\n", "mySearchlight = sl.load('example_data/searchlight_volume.h5')\n", "print(f'Average size of the searchlight: {mySearchlight.radius.mean():.1f} mm')\n", "print(f'Average number of voxels: {mySearchlight.nvoxels.mean():.1f}')" ] }, { "cell_type": "markdown", "id": "7d3e0329", "metadata": {}, "source": [ "To run a searchlight we need:\n", "\n", "* an MVPA function\n", "* datafiles\n", "\n", "The searchlight will calculate the MVPA function on data from datafiles in each center defined in `mySearchlight`.\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "32dbd891", "metadata": {}, "outputs": [], "source": [ "# Here we use simple example that calculates the mean activity in the searchlight region.\n", "# The MVPA functions take a `n_images x n_voxels` array as input and return either a scalar or vector as output\n", "\n", "def mvpa_mean_function(data):\n", " # Example of MVPA function that returns a scalar as an output argument\n", " return np.mean(data)\n", "\n", "def mvpa_multi_function(data):\n", " # Example of MVPA function that returns a vector as an output argument\n", " return np.mean(data,axis=1)\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "b48fa723", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Run the searchlight with the defined scalar MVPA function\n", "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)]\n", "results = mySearchlight.run(datafiles, mvpa_mean_function,verbose=False)\n", "mySearchlight.data_to_nifti(results,outfilename='example_data/output_scalar.nii')" ] }, { "cell_type": "code", "execution_count": 10, "id": "ccdafd87", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Run the searchlight with the defined scalar vector function\n", "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)]\n", "results = mySearchlight.run(datafiles, mvpa_multi_function,verbose=False)\n", "mySearchlight.data_to_nifti(results,outfilename='example_data/output_vector.nii')" ] }, { "cell_type": "markdown", "id": "1ed5f2e1", "metadata": {}, "source": [ "Running representational similarity analysis (RSA) as the MVPA function:" ] }, { "cell_type": "markdown", "id": "e604cada", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 11, "id": "ac0d3c14", "metadata": {}, "outputs": [], "source": [ "# append resms to datafiles\n", "datafiles.append('example_data/sub-02_ses-s1_resms.nii')\n", "\n", "# define function to calculate\n", "def mvpa_mean_rsa(data,cond_vector, partition_vector):\n", " # 1.prewhite the betas\n", " y = np.zeros_like(data)\n", " for i in range(data.shape[0] - 1):\n", " resms = data[-1,:]\n", " resms[resms <= 0] = np.nan\n", " y[i,:] = data[i,:] / np.sqrt(resms)\n", " y = y[:-1,:] # remove the resms from the data\n", " \n", " # 2.Calculate the RDM using rsatoolbox\n", " 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})\n", " rdm = rsa.rdm.calc_rdm(dataset, method='crossnobis', descriptor='conds', cv_descriptor='run')\n", " rdm.rdm_descriptors = {'index': [0]}\n", " \n", " # 3.get the upper triangular\n", " rdm_vec = rdm.dissimilarities\n", "\n", " return np.mean(rdm_vec)\n" ] }, { "cell_type": "markdown", "id": "d4e19cbd", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 12, "id": "dc38a429", "metadata": {}, "outputs": [], "source": [ "condition_vector = [1,2,1,2] # condition number for each beta\n", "partition_vector = [1,1,2,2] # run number for each beta" ] }, { "cell_type": "code", "execution_count": 13, "id": "0945273d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Loading data \n", "Volume number 1\n", "Volume number 2\n", "Volume number 3\n", "Volume number 4\n", "Volume number 5\n", "Calculating searchlight 0 of 70752\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Calculating searchlight 1000 of 70752\n", "Calculating searchlight 2000 of 70752\n", "Calculating searchlight 3000 of 70752\n", "Calculating searchlight 4000 of 70752\n", "Calculating searchlight 5000 of 70752\n", "Calculating searchlight 6000 of 70752\n", "Calculating searchlight 7000 of 70752\n", "Calculating searchlight 8000 of 70752\n", "Calculating searchlight 9000 of 70752\n", "Calculating searchlight 10000 of 70752\n", "Calculating searchlight 11000 of 70752\n", "Calculating searchlight 12000 of 70752\n", "Calculating searchlight 13000 of 70752\n", "Calculating searchlight 14000 of 70752\n", "Calculating searchlight 15000 of 70752\n", "Calculating searchlight 16000 of 70752\n", "Calculating searchlight 17000 of 70752\n", "Calculating searchlight 18000 of 70752\n", "Calculating searchlight 19000 of 70752\n", "Calculating searchlight 20000 of 70752\n", "Calculating searchlight 21000 of 70752\n", "Calculating searchlight 22000 of 70752\n", "Calculating searchlight 23000 of 70752\n", "Calculating searchlight 24000 of 70752\n", "Calculating searchlight 25000 of 70752\n", "Calculating searchlight 26000 of 70752\n", "Calculating searchlight 27000 of 70752\n", "Calculating searchlight 28000 of 70752\n", "Calculating searchlight 29000 of 70752\n", "Calculating searchlight 30000 of 70752\n", "Calculating searchlight 31000 of 70752\n", "Calculating searchlight 32000 of 70752\n", "Calculating searchlight 33000 of 70752\n", "Calculating searchlight 34000 of 70752\n", "Calculating searchlight 35000 of 70752\n", "Calculating searchlight 36000 of 70752\n", "Calculating searchlight 37000 of 70752\n", "Calculating searchlight 38000 of 70752\n", "Calculating searchlight 39000 of 70752\n", "Calculating searchlight 40000 of 70752\n", "Calculating searchlight 41000 of 70752\n", "Calculating searchlight 42000 of 70752\n", "Calculating searchlight 43000 of 70752\n", "Calculating searchlight 44000 of 70752\n", "Calculating searchlight 45000 of 70752\n", "Calculating searchlight 46000 of 70752\n", "Calculating searchlight 47000 of 70752\n", "Calculating searchlight 48000 of 70752\n", "Calculating searchlight 49000 of 70752\n", "Calculating searchlight 50000 of 70752\n", "Calculating searchlight 51000 of 70752\n", "Calculating searchlight 52000 of 70752\n", "Calculating searchlight 53000 of 70752\n", "Calculating searchlight 54000 of 70752\n", "Calculating searchlight 55000 of 70752\n", "Calculating searchlight 56000 of 70752\n", "Calculating searchlight 57000 of 70752\n", "Calculating searchlight 58000 of 70752\n", "Calculating searchlight 59000 of 70752\n", "Calculating searchlight 60000 of 70752\n", "Calculating searchlight 61000 of 70752\n", "Calculating searchlight 62000 of 70752\n", "Calculating searchlight 63000 of 70752\n", "Calculating searchlight 64000 of 70752\n", "Calculating searchlight 65000 of 70752\n", "Calculating searchlight 66000 of 70752\n", "Calculating searchlight 67000 of 70752\n", "Calculating searchlight 68000 of 70752\n", "Calculating searchlight 69000 of 70752\n", "Calculating searchlight 70000 of 70752\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result = mySearchlight.run(datafiles,mvpa_mean_rsa,function_args={'cond_vector': condition_vector, 'partition_vector': partition_vector},verbose=True)\n", "mySearchlight.data_to_nifti(result,outfilename='example_data/output_scalar_RSA.nii')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.0" } }, "nbformat": 4, "nbformat_minor": 5 }