Cortical Hemisphere

Here, we’ll perform an analysis similar to those included in our pre-print, using our example data.

First, create an instance of the Base class and generate 1000 surrogate maps:

from brainsmash.mapgen.base import Base
import numpy as np

# load parcellated structural neuroimaging maps
myelin = np.loadtxt("LeftParcelMyelin.txt")
thickness = np.loadtxt("LeftParcelThickness.txt")

# instantiate class and generate 1000 surrogates
gen = Base(myelin, "LeftParcelGeodesicDistmat.txt")  # note: can pass numpy arrays as well as filenames
surrogate_maps = gen(n=1000)

Next, compute the Pearson correlation between each surrogate myelin map and the empirical cortical thickness map:

from brainsmash.mapgen.stats import pearsonr, pairwise_r

surrogate_brainmap_corrs = pearsonr(thickness, surrogate_maps).flatten()
surrogate_pairwise_corrs = pairwise_r(surrogate_maps, flatten=True)

Repeat using randomly shuffled surrogate myelin maps:

naive_surrogates = np.array([np.random.permutation(myelin) for _ in range(1000)])
naive_brainmap_corrs = pearsonr(thickness, naive_surrogates).flatten()
naive_pairwise_corrs = pairwise_r(naive_surrogates, flatten=True)

Now plot the results:

import matplotlib.pyplot as plt
from scipy import stats

sac = '#377eb8'  # autocorr-preserving
rc = '#e41a1c'  # randomly shuffled
bins = np.linspace(-1, 1, 51)  # correlation b

# this is the empirical statistic we're creating a null distribution for
test_stat = stats.pearsonr(myelin, thickness)[0]

fig = plt.figure(figsize=(3, 3))
ax = fig.add_axes([0.2, 0.25, 0.6, 0.6])  # autocorr preserving
ax2 = ax.twinx()  # randomly shuffled

# plot the data
ax.axvline(test_stat, 0, 0.8, color='k', linestyle='dashed', lw=1)
ax.hist(surrogate_brainmap_corrs, bins=bins, color=sac, alpha=1,
    density=True, clip_on=False, zorder=1)
ax2.hist(naive_brainmap_corrs, bins=bins, color=rc, alpha=0.7,
    density=True, clip_on=False, zorder=2)

# make the plot nice...
ax.set_xticks(np.arange(-1, 1.1, 0.5))
ax.tick_params(axis='y', colors=sac)
ax2.tick_params(axis='y', colors=rc)
ax.set_ylim(0, 2)
ax2.set_ylim(0, 6)
ax.set_xlim(-1, 1)
[s.set_visible(False) for s in [
    ax.spines['top'], ax.spines['right'], ax2.spines['top'], ax2.spines['left']]]
ax.text(0.97, 1.1, 'SA-independent', ha='right',va='bottom',
    color=rc, transform=ax.transAxes)
ax.text(0.97, 1.03, 'SA-preserving', ha='right', va='bottom',
    color=sac, transform=ax.transAxes)
ax.text(test_stat, 1.65, "T1w/T2w\nmap", ha='center', va='bottom')
ax.text(0.5, -0.2, "Pearson correlation\nwith thickness map",
    ha='center', va='top', transform=ax.transAxes)
ax.text(-0.3, 0.5, "Density", rotation=90, ha='left', va='center', transform=ax.transAxes)

Executing the above code produces the following figure:


We can plot a couple surrogate maps on the cortical surface using wbplot:

from wbplot import pscalar

def vrange(x):
    return (np.percentile(x, 5), np.percentile(x, 95))

for i in range(3):
    y = surrogate_maps[i]

Executing the above code produces the following three images:


We’ll assess our surrogate maps’ reliability using their fit to the parcellated myelin map’s variogram:

from brainsmash.mapgen.eval import base_fit

    nh=25,  # these are default kwargs, but shown here for demonstration
    deltas=np.arange(0.1, 1, 0.1),
    pv=25)  # kwargs are passed to brainsmash.mapgen.base.Base

Executing the code above produces the following plot:


The surrogate maps exhibit the same autocorrelation structure as the empirical brain map.

Finally, we’ll compute non-parametric p-values using our two different null distributions:

from brainsmash.mapgen.stats import nonparp

print("Spatially naive p-value:", nonparp(test_stat, naive_brainmap_corrs))
print("SA-corrected p-value:", nonparp(test_stat, surrogate_brainmap_corrs))

The two p-values for this example come out to P < 0.001 and P=0.001, respectively.

Unilateral Subcortical Structure

For a subcortical analysis you’ll typically need:

  • A subcortical distance matrix

  • A subcortical brain map

  • A mask corresponding to a structure of interest

We’ll assume you use Connectome Workbench-style files and that you want to isolate one particular anatomical structure.


If you already have a distance matrix and a brain map for your subcortical structure of interest, the workflow is identical to the cortical examples in Getting Started.

If you haven’t already computed a subcortical distance matrix or downloaded our pre-computed version, please follow these steps.

To isolate one subcortical structure from a whole-brain dscalar file, first do:

wb_command -cifti-export-dense-mapping yourfile.dscalar.nii COLUMN -volume-all output.txt -structure

We can then use the information contained in this file to isolate one particular structure, e.g. the left cerebellum:

from brainsmash.utils.dataio import load
import numpy as np
import pandas as pd

# Input files
image = "/path/to/yourfile.dscalar.nii"
wb_output = "output.txt"

# Load the output of the above command
df = pd.read_table(wb_output, header=None, index_col=0, sep=' ',
         names=['structure', 'mni_i', 'mni_j', 'mni_k']).rename_axis('index')

# Get indices for left cerebellum
indices = df[df['structure'] == 'CEREBELLUM_LEFT'].index.values

# Create a binary mask
mask = np.ones(31870)  # volume has 31870 CIFTI indices in standard 91k mesh
indices -= 59412  # first 59412 indices in whole-brain maps are cortical
mask[indices] = 0
np.savetxt("mask.txt", mask)  # this mask has right dimension for distmat

# Also saved a masked copy of the image
image_data = load(image)
indices += 59412  # assuming image data is whole-brain!
masked_image = image_data[indices]
np.savetxt("masked_image.txt", masked_image)

Next, we’ll need to sort and memory-map our distance matrix, but only for the pairwise distances between left cerebellar voxels:

from brainsmash.mapgen.memmap import txt2memmap

# Input files
image = "masked_image.txt"
mask = "mask.txt"
distmat = "/path/to/SubcortexDenseEuclideanDistmat.txt"

output_files = txt2memmap(distmat, output_dir=".", maskfile=mask, delimiter=' ')

Now, we can use the output files to instantiate our surrogate map generator. Here, we’ll also use the keyword arguments which were used in our study for left cerebellum. First, we’ll validate the variogram fit using these parameters:

from brainsmash.mapgen.eval import sampled_fit

x = "masked_image.txt"
distmat = output_files['distmat']
index = output_files['index']

kwargs = {'ns': 500,
          'knn': 1500,
          'nh': 25,
          'deltas': [0.3, 0.5, 0.7, 0.9],
          'pv': 70

sampled_fit(x, distmat, index, nsurr=20, **kwargs)

This produces the following plot:


Having confirmed that the fit looks good, we simulate cerebellar surrogate maps with a call to the surrogate map generator:

from brainsmash.mapgen.sampled import Sampled

gen = Sampled(x, distmat, index, **kwargs)
surrogate_maps = gen(n=100)

Whole-Brain Volume

For this analysis you will need:

  • A text file containing voxel coordinates (with shape N rows by 3 columns)

  • A text file containing N brain map values

Your favorite neuroimaging software should be able to generate these files (e.g., AFNI).


Do not include every voxel in the volume, but rather the voxels corresponding to your region of interest (e.g., the brain).

First, you will need to generate memory-mapped distance matrix files. This is achieved using brainsmash.workbench.geo.volume, as in the code below:

from brainsmash.workbench.geo import volume

coord_file = "/path/to/voxel_coordinates.txt"
output_dir = "/some/directory/"

filenames = volume(coord_file, output_dir)

A dictionary containing the names of the newly generated files is returned by the function, which are used in the example below. Prior to generating surrogate maps, we first visually inspect the variogram fit:

from brainsmash.mapgen.eval import sampled_fit

brain_map = "/some_path_to/brain_map.txt"

# These are three of the key parameters affecting the variogram fit
kwargs = {'ns': 500,
          'knn': 1500,
          'pv': 70

# Running this command will generate a matplotlib figure
sampled_fit(brain_map, filenames['D'], filenames['index'], nsurr=10, **kwargs)

If you are unhappy with the fit, it is recommended that you try varying the three parameters included in the keyword argument dictionary above (i.e., ns, knn, and pv).


A few users have reported a systematic overestimation of variance in their surrogate maps’ variograms. This suggests that further parameter tuning is necessary. However, it may also indicate a violation of the underlying assumption of stationary, which is more likely to be violated at the whole-brain level than within a single contiguous brain structure.

Having selected our parameter values and stored them in the kwargs dictionary, we can now randomly generate surrogate brain maps with spatial autocorrelation that is matched to our target brain map:

from brainsmash.mapgen.sampled import Sampled

gen = Sampled(x=brain_map, D=filenames['D'], index=filenames['index'], **kwargs)
surrogate_maps = gen(n=1000)

These surrogate maps can then be used to test the null hypothesis that any random yet spatially autocorrelated brain map would yield a comparable or more extreme statistical result.