Building gradients

This section describes how to generate macroscale gradient mapping from the output matrices of micapipe. The matrices are the same as in the Main output matrices tutorial.

For this tutorial we will map each modality of a single subject using BrainSpace, a python based library. Additionally the libraries os, numpy, nibabel and nibabel will be used.

As in previous examples, we will use the subject HC001, session 01 from the MICs dataset, and all paths will be relative to the subject directory or out/micapipe/sub-HC001_ses01/ and the atlas schaefer-400.

python environment

The first step is to set the python environment and all the variables relative to the subject you’re interested to visualize.

 1# Set the environment
 2import os
 3import numpy as np
 4import nibabel as nb
 5from brainspace.plotting import plot_hemispheres
 6from brainspace.mesh.mesh_io import read_surface
 7from brainspace.datasets import load_conte69
 8from brainspace.gradient import GradientMaps
 9from brainspace.utils.parcellation import map_to_labels
10import matplotlib.pyplot as plt
11
12# Set the working directory to the 'out' directory
13os.chdir("~/out")     # <<<<<<<<<<<< CHANGE THIS PATH
14
15# Path to micapipe repository
16micapipe='~/micapipe' # <<<<<<<<<<<< CHANGE THIS PATH
17
18# This variable will be different for each subject
19subjectID='sub-HC001_ses-01'           # <<<<<<<<<<<< CHANGE THIS SUBJECT's ID
20subjectDir='micapipe/sub-HC001/ses-01' # <<<<<<<<<<<< CHANGE THIS SUBJECT's DIRECTORY
21
22# Here we define the atlas
23atlas='schaefer-400' # <<<<<<<<<<<< CHANGE THIS ATLAS

Loading the surfaces

 1 # Load conte69
 2 c69_lh, c69_rh = load_conte69()
 3
 4 # Load fsaverage5
 5 fs5_lh = read_surface('freesurfer/fsaverage5/surf/lh.pial', itype='fs')
 6 fs5_rh = read_surface('freesurfer/fsaverage5/surf/rh.pial', itype='fs')
 7
 8 # Load annotation file in fsaverage5
 9 annot_lh_fs5= nb.freesurfer.read_annot(micapipe + '/parcellations/lh.schaefer-400_mics.annot')
10 annot_rh_fs5= nb.freesurfer.read_annot(micapipe + '/parcellations/rh.schaefer-400_mics.annot')[0]+200
11
12 # replace with 0 the medial wall of the right labels
13 annot_rh_fs5 = np.where(annot_rh_fs5==200, 0, annot_rh_fs5)
14
15 # fsaverage5 labels
16 labels_fs5 = np.concatenate((annot_lh_fs5[0], annot_rh_fs5), axis=0)
17
18 # Read label for conte69
19 labels_c69 = np.loadtxt(open(micapipe + '/parcellations/schaefer-400_conte69.csv'), dtype=np.int)

Functional gradients

Load and slice the functional matrix

 1 # Load the cortical connectome
 2 mtx_fs = np.loadtxt(subjectDir + '/func/surfaces/' + subjectID + '_rsfmri_space-fsnative_atlas-' + atlas + '_desc-FC.txt',
 3                     dtype=np.float, delimiter=' ')
 4
 5 # slice the matrix
 6 FC = mtx_fs[49:, 49:]
 7 FC = np.delete(np.delete(FC, 200, axis=0), 200, axis=1)
 8
 9 # Fischer transformation
10 FCz = np.arctanh(FC)
11
12 # replace inf with 0
13 FCz[~np.isfinite(FCz)] = 0
14
15 # Mirror the matrix
16 FCz = np.triu(FCz,1)+FCz.T

Calculate the functional gradients

1 # Number of gradients
2 N = 10
3
4 # Calculate the gradients
5 gm = GradientMaps(n_components=N, random_state=None, approach='dm', kernel='normalized_angle')
6 gm.fit(FCz, sparsity=0.8)

Plot the functional gradients

 1 # Plot the gradients
 2 g1=gm.gradients_[:, 0]
 3 g2=gm.gradients_[:, 1]
 4 g3=gm.gradients_[:, 2]
 5
 6 # Creating figure
 7 fig = plt.subplots(1, 2, figsize = (7, 5))
 8 ax = plt.axes(projection ="3d")
 9
10 # Creating plot
11 ax.scatter3D(g1, g2, g3, color = 'red')
12 plt.title("Functional gradient")
13 ax.set_xlabel('Grad 1')
14 ax.set_ylabel('Grad 2')
15 ax.set_zlabel('Grad 3')
16
17 # show plot
18 plt.show()
../../_images/fc_scatter.png

Functional gradients to fsaverage5 surface

 1 # Mask of the medial wall on fsaverage 5
 2 mask_fs5 = labels_fs5 != 0
 3
 4 # Map gradients to original parcels
 5 grad = [None] * 3
 6 for i, g in enumerate(gm.gradients_.T[0:3,:]):
 7     grad[i] = map_to_labels(g, labels_fs5,  fill=np.nan, mask=mask_fs5)
 8
 9 # Plot Gradients RdYlBu
10 plot_hemispheres(fs5_lh, fs5_rh, array_name=grad, size=(1000, 600), cmap='coolwarm',
11                  embed_nb=True,  label_text={'left':['Grad1', 'Grad2','Grad3']}, color_bar='left',
12                  zoom=1.25, nan_color=(1, 1, 1, 1) )
../../_images/fc_fs5.png

Functional gradients to conte69 surface

 1 # mask of the medial wall
 2 mask_c69 = labels_c69 != 0
 3
 4 # Map gradients to original parcels
 5 grad = [None] * 3
 6 for i, g in enumerate(gm.gradients_.T[0:3,:]):
 7     grad[i] = map_to_labels(g, labels_c69, fill=np.nan, mask=mask_c69)
 8
 9 # Plot Gradients coolwarm
10 plot_hemispheres(c69_lh, c69_rh, array_name=grad, size=(1000, 600), cmap='coolwarm',
11                  embed_nb=True,  label_text={'left':['Grad1', 'Grad2','Grad3']}, color_bar='left',
12                  zoom=1.25, nan_color=(1, 1, 1, 1))
../../_images/fc_c69.png

Structural gradients

Load and slice the structural matrix

 1 # Load the cortical connectome
 2 mtx_sc = np.loadtxt(subjectDir + '/dwi/connectomes/' + subjectID + '_space-dwi_atlas-' + atlas + '_desc-iFOD2-40M-SIFT2_cor-connectome.txt',
 3                     dtype=np.float, delimiter=' ')
 4
 5 # Fill the lower triangle of the matrix
 6 mtx_sc = np.log(np.triu(mtx_sc,1)+mtx_sc.T)
 7 mtx_sc[np.isneginf(mtx_sc)] = 0
 8
 9 # Slice the connectome to use only cortical nodes
10 SC = mtx_sc[49:, 49:]
11 SC = np.delete(np.delete(SC, 200, axis=0), 200, axis=1)

Calculate the structural gradients

1 # SC Left hemi
2 gm_SC_L = GradientMaps(n_components=N, random_state=None, approach='dm', kernel='normalized_angle')
3 gm_SC_L.fit(SC[0:200, 0:200], sparsity=0)
4
5 # SC Right hemi
6 gm_SC_R = GradientMaps(n_components=N, alignment='procrustes', kernel='normalized_angle'); # align right hemi to left hemi
7 gm_SC_R.fit(SC[200:400, 200:400], sparsity=0, reference=gm_SC_L.gradients_)

Plot the structural gradients

 1 # plot the left gradients
 2 g1=gm_SC_L.gradients_[:, 0]
 3 g2=gm_SC_L.gradients_[:, 1]
 4 g3=gm_SC_L.gradients_[:, 2]
 5 # plot the right gradients
 6 g1R=gm_SC_R.gradients_[:, 0]
 7 g2R=gm_SC_R.gradients_[:, 1]
 8 g3R=gm_SC_R.gradients_[:, 2]
 9
10 # Creating figure
11 fig = plt.subplots(1, 2, figsize = (7, 5))
12 ax = plt.axes(projection ="3d")
13
14 # Creating plot
15 ax.scatter3D(g1, g2, g3, color = 'purple')
16 ax.scatter3D(g1R, g2R, g3R, color = 'slateblue', marker='v')
17 plt.title("Structural gradient")
18 ax.legend(['Left SC', 'Right SC'])
19 ax.set_xlabel('Grad 1')
20 ax.set_ylabel('Grad 2')
21 ax.set_zlabel('Grad 3')
22
23 # show plot
24 plt.show()
../../_images/sc_scatter.png

Structural gradients to conte69 surface

 1 # Left and right gradients concatenated
 2 SC_gradients = np.concatenate((gm_SC_L.gradients_, gm_SC_R.gradients_), axis=0)
 3
 4 # Map gradients to original parcels
 5 grad = [None] * 3
 6 for i, g in enumerate(SC_gradients.T[0:3,:]):
 7     grad[i] = map_to_labels(g, labels_c69, fill=np.nan, mask=mask_c69)
 8
 9 # Plot Gradients
10 plot_hemispheres(c69_lh, c69_rh, array_name=grad, size=(1000, 600), cmap='coolwarm',
11                  embed_nb=True,  label_text={'left':['Grad1', 'Grad2','Grad3']}, color_bar='left',
12                  zoom=1.25, nan_color=(1, 1, 1, 1) )
../../_images/sc_c69.png

MPC gradients

Load and slice the MPC matrix

 1 # Load the cortical connectome
 2 mtx_mpc = np.loadtxt(subjectDir + '/anat/surfaces/micro_profiles/' + subjectID + '_space-fsnative_atlas-' + atlas + '_desc-MPC.txt',
 3                      dtype=np.float, delimiter=' ')
 4
 5 # Fill the lower triangle of the matrix
 6 MPC = np.triu(mtx_mpc,1)+mtx_mpc.T
 7
 8 # Renove the medial wall
 9 MPC = np.delete(np.delete(MPC, 0, axis=0), 0, axis=1)
10 MPC = np.delete(np.delete(MPC, 200, axis=0), 200, axis=1)

Calculate the MPC gradients

1 # Calculate the gradients
2 gm = GradientMaps(n_components=15, random_state=None, approach='dm', kernel='normalized_angle')
3 gm.fit(MPC, sparsity=0.8)

Plot the MPC gradients

 1 # Plot the gradients
 2 g1=gm.gradients_[:, 0]
 3 g2=gm.gradients_[:, 1]
 4 g3=gm.gradients_[:, 2]
 5
 6 # Creating figure
 7 fig = plt.subplots(1, 2, figsize = (7, 5))
 8 ax = plt.axes(projection ="3d")
 9
10 # Creating plot
11 ax.scatter3D(g1, g2, g3, color = 'green')
12 plt.title("MPC gradient")
13 ax.set_xlabel('Grad 1')
14 ax.set_ylabel('Grad 2')
15 ax.set_zlabel('Grad 3')
16
17 # show plot
18 plt.show()
../../_images/mpc_scatter.png

MPC gradients to conte69 surface

1 # Map gradients to original parcels
2 grad = [None] * 3
3 for i, g in enumerate(gm.gradients_.T[0:3,:]):
4     grad[i] = map_to_labels(g, labels_c69, fill=np.nan, mask=mask_c69)
5
6 # Plot Gradients
7 plot_hemispheres(c69_lh, c69_rh, array_name=grad, size=(1000, 600), cmap='coolwarm',
8                  embed_nb=True,  label_text={'left':['MPC-G1', 'MPC-G2','MPC-G3']}, color_bar='left',
9                  zoom=1.25, nan_color=(1, 1, 1, 1) )
../../_images/mpc_c69.png

Geodesic distance gradients

Load and slice the GD matrix

1 # Load the cortical connectome
2 mtx_gd = np.loadtxt(subjectDir + '/anat/surfaces/geo_dist/' + subjectID + '_space-fsnative_atlas-' + atlas + '_GD.txt',
3                     dtype=np.float, delimiter=' ')
4
5 # Remove the Mediall Wall
6 mtx_gd = np.delete(np.delete(mtx_gd, 0, axis=0), 0, axis=1)
7 GD = np.delete(np.delete(mtx_gd, 200, axis=0), 200, axis=1)

Calculate the GD gradients

1 # GD Left hemi
2 gm_GD_L = GradientMaps(n_components=N, random_state=None, approach='dm', kernel='normalized_angle')
3 gm_GD_L.fit(GD[0:200, 0:200], sparsity=0.8)
4
5 # GD Right hemi
6 gm_GD_R = GradientMaps(n_components=N, alignment='procrustes', kernel='normalized_angle'); # align right hemi to left hemi
7 gm_GD_R.fit(GD[200:400, 200:400], sparsity=0.8, reference=gm_GD_L.gradients_)

Plot the GD gradients

 1 # plot the gradients
 2 g1=gm_GD_L.gradients_[:, 0]
 3 g2=gm_GD_L.gradients_[:, 1]
 4 g3=gm_GD_L.gradients_[:, 2]
 5 # plot the gradients
 6 g1R=gm_GD_R.gradients_[:, 0]
 7 g2R=gm_GD_R.gradients_[:, 1]
 8 g3R=gm_GD_R.gradients_[:, 2]
 9
10 # Creating figure
11 fig = plt.subplots(1, 2, figsize = (7, 5))
12 ax = plt.axes(projection ="3d")
13
14 # Creating plot
15 ax.scatter3D(g1, g2, g3, color = 'dodgerblue')
16 ax.scatter3D(g1R, g2R, g3R, color = 'teal', marker='v')
17 plt.title("Structural gradient")
18 ax.legend(['Left GD', 'Right GD'])
19 ax.set_xlabel('Grad 1')
20 ax.set_ylabel('Grad 2')
21 ax.set_zlabel('Grad 3')
22
23 # show plot
24 plt.show()
../../_images/gd_scatter.png

GD gradients to conte69 surface

 1 # Left and right gradients concatenated
 2 GD_gradients = np.concatenate((gm_GD_L.gradients_, gm_GD_R.gradients_), axis=0)
 3
 4 # Map gradients to original parcels
 5 grad = [None] * 3
 6 for i, g in enumerate(GD_gradients.T[0:3,:]):
 7     grad[i] = map_to_labels(g, labels_c69, fill=np.nan, mask=mask_c69)
 8
 9 # Plot Gradients
10 plot_hemispheres(c69_lh, c69_rh, array_name=grad, size=(1000, 600), cmap='coolwarm',
11                  embed_nb=True,  label_text={'left':['GD-G1', 'GD-G1','GD-G3']}, color_bar='left',
12                  zoom=1.25, nan_color=(1, 1, 1, 1))
../../_images/gd_c69.png

Download the code!

Python Jupyter notebook: 'tutorial_gradients.ipynb'

Python source code: 'tutorial_gradients.py'