Spin test
The spin test is a permutation technique designed to assess the spatial correspondence between two brain maps of neuroanatomical structure-specific surface measures, whilst accounting for the possibility of two surfaces being randomly aligned.
It was developed by Alexander-Bloch et al., 2018, and is compared to other techniques for assessing spatial relationships between brain maps in Markello et al., 2020.
The test works by generating a null distribution of a correspondence statistic between randomly permuted versions ("spins") of the first brain map and the original version of the second brain map. The output p-value corresponds to the proportion of randomly generated correspondence estimates that are above the observed correspondence between the two brain maps.
Thus, one first computes the observed correspondence statistic (below, we were looking at the Pearson correlation coefficient) between the two brain maps of interest. Then, the first brain map is projected onto a sphere, spatially permuted by a random angle, and projected back onto the structure of interest. The correspondence statistic is then computed between the spatially permutated brain map and the second, untouched brain map. This process is repeated n times, in order to create a distribution of estimates of the correspondence statistic.
First, one needs to decide on a correspondence statistic to examine between the two brain maps. This choice depends on the data and scientific question of interest.
For instance, the normalized mutual information (NMI) coefficient can be selected in the case of two brains maps of categorical data (e.g. comparing cortical parcellations, like in Alexander-Bloch et al., 2018).
For the following wiki page, we will be working under the assumption that we are comparing two brain maps of continuous data points and have therefore selected the Pearson's correlation coefficient as the correspondence statistic.
Modifications to the original spin test scripts were performed in Matlab R2020b (9.9.0). Matlab is available for free through McGill and requires 10 GB of free space.
It is not available on the CIC, so you will need to work from your local computer.
The first thing to do is to clone the original spin-test GitHub repository:
git clone https://github.com/spin-test/spin-test.git
Here are the 3 most important scripts in the spin-test-master directory:
- SpinPermuCIVET.m : projects the vertex files onto a sphere object and spins one surface the specified number of times. Outputs the spinned surfaces
- pvalvsNull.m : generates the spatial correspondence statistics (e.g. correlation coefficient) for the original and spinned surfaces and compares the statistics. Outputs the original correspondence statistic and p-value
- run_spintest.m : specifies the file paths (inputs and outputs) and runs the two previous scripts
You will be replacing the pvalvsNull.m and run_spintest.m scripts with the modified version to be found on the following wiki:
NOTE: This section is currently under construction and should not be used. The order of vertices along a custom sphere generated using this method does not recapitulate the order of data in surface (vertex) text files, e.g. from CIVET.
Input: surface_of_interest.obj (ipsilateral surface of interest)
Output: sphere.obj (spherical object with same number of vertices as surface_of_interest.obj)
Note: script below is a collage from StackOverflow.
head -1 surface_of_interest.obj # the number of vertices will be the last integer
Make sure to replace N on line 5 with the number of vertices you found in step 1
# PYTHON 3.8.8
! pip install vtk # optional
import numpy as np
import vtk
num_pts = N
indices = np.arange(0, num_pts, dtype=float) + 0.5
phi = np.arccos(1 - 2*indices/num_pts)
theta = np.pi * (1 + 5**0.5) * indices
x, y, z = np.cos(theta) * np.sin(phi), np.sin(theta) * np.sin(phi), np.cos(phi);
points = vtk.vtkPoints()
for i in range(len(x)):
array_point = np.array([x[i], y[i], z[i]] )
points.InsertNextPoint(x[i],y[i],z[i])
poly = vtk.vtkPolyData()
poly.SetPoints(points)
# To create surface of a sphere we need to use Delaunay triangulation
d3D = vtk.vtkDelaunay3D()
d3D.SetInputData( poly ) # This generates a 3D mesh
# We need to extract the surface from the 3D mesh
dss = vtk.vtkDataSetSurfaceFilter()
dss.SetInputConnection( d3D.GetOutputPort() )
dss.Update()
# Now we have our final polydata
spherePoly = dss.GetOutput()
writer = vtk.vtkPolyDataWriter()
writer.SetInputData(spherePoly)
writer.SetFileName("sphere.vtk")
writer.Write()
module load minc-toolkit/1.9.11
module load anaconda/5.1.0-python3
module load pyminc/0.51
module load minc-stuffs/0.1.9^minc-toolkit-1.9.10
vtk_meshconvert.py -i sphere.vtk -a -o sphere.obj
Replace the sphere.obj file in your spin-test-master/SurfStat directory with the one that you just created.
Added a linear regression on the z-scored surface maps.
-
This then allows for the visualization of the relationship in graph form and the calculation of residuals of that relationship.
-
The residuals can then be mapped back to the cortical surface and visualized to assess which areas exhibit the observed relationship (i.e. lower residuals), and which areas do not (i.e. higher residuals).
NOTE: If you are using non-cortical data, you can comment out or modify lines 10-16 (designed to exclude the cortical medial well), depending on whether or not you have vertices to exclude and what those vertices are.
If you are working with cortical data, it is highly recommended to exclude the medial wall. To do so, you can download the left and right medial wall masks from the CIC (/data/chamal/projects/paroli/ressources/CIVET_2.0_mask_left_short_inv.txt and /data/chamal/projects/paroli/ressources/CIVET_2.0_mask_right_short_inv.txt), and move these files into the "data" folder of your spin test directory
function pval=pvalvsNull_resid(x_label,y_label,readleft1,readright1,readleft2,readright2,permno,wsname, real_rho_path, resid_path_left, resid_path_right)
%load the saved workspace from SpinPermu
load(wsname);
%read the data saved in csv and merge left and right surfaces into one
datal1=importdata(readleft1);
datar1=importdata(readright1);
datal2=importdata(readleft2);
datar2=importdata(readright2);
% Label medial wall vertices with NaN:
v_exclude_left = logical(load('../data/CIVET_2.0_mask_left_short_inv.txt')); % change path
v_exclude_right = logical(load('../data/CIVET_2.0_mask_right_short_inv.txt'));
datal1(v_exclude_left) = NaN;
datal2(v_exclude_left) = NaN;
datar1(v_exclude_right) = NaN;
datar2(v_exclude_right) = NaN;
data1=cat(1,datal1,datar1);
data2=cat(1,datal2,datar2);
% Z-score input data
zscor_xnan = @(x) bsxfun(@rdivide, bsxfun(@minus, x, mean(x,'omitnan')), std(x, 'omitnan')); % Function to Z-score surfaces with the NaNs from medial wall removal
Z_data1 = zscor_xnan(data1); %Z-scored values of first map
Z_data2 = zscor_xnan(data2); %Z-scored values of second map
%calculate the real Pearson's correlation between two interested maps
realrho=corr(data1,data2, 'rows','complete'); % 'rows','complete' to exclude NaN's
Z_realrho = corr(Z_data1,Z_data2, 'rows','complete'); %Sanity check to see if correlation with original vs z-scored data is the same (it should be)
% Get residuals from regression
[b,bint,resid] = regress(Z_data1,Z_data2); % Get residuals from regression of z-scored surfaces
resid_left = resid(1:40962,:); % Divide into left and right
resid_right = resid(40963:81924,:);
resid_left(v_exclude_left) = 0; %residual values on medial wall = 0
resid_right(v_exclude_right) = 0;
dlmwrite(strcat(output_dir,'left_residuals.txt'), resid_left); %write residuals as vertex files (.txt)
dlmwrite(strcat(output_dir,'right_residuals.txt'), resid_right);
% Scatter plot of vertices and regression line
model = fitlm(Z_data1,Z_data2); % Regression function on z-scored surfaces
beta_graph_reg = table2array(model.Coefficients(2,1)); %graph regression + scatter plot
figure('visible','off') %so the graph doesn't pop up when running the code
plot(model); %generate graph
graph_path = strcat(output_dir, '/vertex_residual_plot.png'); %where to save graph
saveas(gcf, graph_path) %save graph as .png
% test the observed rho against null described by SpinPermu
nullrho=[];
for i=1:permno
tempdata=cat(2,bigrotl(i,:),bigrotr(i,:))';
nullrho=cat(1,nullrho,corr(tempdata,data2, 'rows','complete')); % 'rows','complete' to exclude NaN's
end
%assuming sign is preserved, calculate the probability that the observed
%correlation coeffcient is above the null distribution
pval=length(find(abs(nullrho)>abs(realrho)))/permno; % added abs() 07/31/2020pval
writematrix(nullrho,sprintf('%s/PearsonRs.csv',output_dir)) % save list of randomly generated correlations
% create & save histogram of randomly generated correlations (can be valuable downstream):
h=figure
hist(nullrho,permno)
xlabel(sprintf("Pearson R"))
title(sprintf('Pearson R= %d; p= %d',realrho,pval))
graph_path = strcat(output_dir,'/correspondence_stat_histogram.png')
saveas(h, graph_path)
dlmwrite(strcat(output_dir,'correspondence_stat_null.txt'), nullrho);
coeffs = [realrho; Z_realrho; b; beta_graph_reg]; %all coefficients should be the same (correlation of original surfaces, correlation of z-scored surfaces, regression on z-scored surfaces for residuals, regression on z-scored surfaces for vertex plot)
dlmwrite(strcat(output_dir,'coefficients.txt'), coeffs);
The scripts are designed to run from the 'spin-test-master' directory
readleft = /path/to/surface/map/to/spin/first_surface_vertex_file_left.csv;
readright = /path/to/surface/map/to/spin/first_surface_vertex_file_right.csv;
permno = 1000; % how many spins (recommended is 1000)
wsname = /path/to/output/spinned/surface/first_surface_rotations.mat;
SpinPermuCIVET(readleft,readright,permno,wsname)
This script will likely run for 20-40 minutes, depending on your computer and the number of permutations specified
readleft1 = readleft;
readright1 = readright;
readleft2 = /path/to/second/surface/map/second_surface_vertex_file_left.csv;
readright2 = /path/to/second/surface/map/second_surface_vertex_file_right.csv;
real_rho_path = /path/to/output/correlation/coefficient/first_second_coef.txt;
resid_path_left = /path/to/output/residuals/first_second_residuals_left.txt;
resid_path_right = /path/to/output/residuals/first_second_residuals_right.txt;
pval=pvalvsNull_resid(readleft1,readright1,readleft2,readright2,permno,wsname,real_rho_path,resid_path_left,resid_path_right);
dlmwrite(/path/to/output/pvalue/first_second_pval.txt), pval);
This script should run in under 1 minute
coefficients.txt
- The correlation coefficients of the original cortical maps;
- There are four values: 1) correlation of original surfaces 2) correlation of z-scored surfaces 3) regression on z-scored surfaces for residuals 4) regression on z-scored surfaces for vertex plot. These values should all be the same (sanity check);
- If value is positive, indicates that where the values are higher on the first surface, values also tend to be higher on the second surface (and vice-versa);
- If value is negative, indicates that where the values are higher on the first surface, values tend to be lower on the second surface (and vice-versa).
correspondence_stat_null.txt
- List of randomly generated estimates for the correspondence statistic;
- You can take a look at the the distribution of that statistic through on the graph saved as
correspondence_stat_histogram.png:
pvalue.txt
- P-value of the correlation coefficient;
- Definition : number of times the correlation coefficient of the spinned surfaces was higher than the original correlation coefficient (in absolute values), divided by the number of permutations;
- This value can be 0, which means that none of the correlations of the spinned surfaces were above the correlation of the original surfaces. In this case, any multiple comparisons correction that divides the p-value by something (e.g. FDR, Bonferonni), will not modify the p-value.
graph.png
- Vertex graph and regression line between the two original surfaces;
- Every point is one vertex, the y-axis is the value on the first surface and the x-axis is the value on the second surface;
- This graph can be used to better interpret the source of the higher residuals;
- It can also be useful to assess the fit of correlation (e.g. if the relationship is more linear, quadratic, cubic...);
residuals_left.txt and residuals_right.txt
- Vertex files of residuals from the regression between the two z-scored surfaces, which can then be visualized with the
create_civet_image.shscript (in minc_extras module on CIC); - Values are expressed in standard deviations;
- Positive values indicate vertices above the regression line, and negative values indicate vertices below the regression line;
- The closer a vertex-wise residual value is to zero, the more that vertex is driving the correlation of interest; it may therefore be interested to visualize the residuals with an inverted colour map.\
Residuals plotted on the cortex with the rminc red and rminc blue colour maps, thresholded at +/- 1 SD. Figure by Olivier Parent.
Residuals plotted on the striatum, thalamus and globus pallidus with an inverted colour map (matplotlib's twilight). Figure by Nadia Blostein.