These are some vary parochial notes solving a particular problem that a colleague of mine had with making surface visualizations of some of their macaque brain data. In other words it may or may not be of general interest...
Notebook Setup¶
Define some variables
# put some system-specific variables in the namespace ('le' dict)
%run ~/set_localenv_vars.py
nb_name = 'macaque_brain_conmat_viz'
# this is where I look for stuf
outdir = '%s/notebooks/%s' %(le['data_dir'],nb_name)
# this is where stuff actually lives
outdir_loc = '%s/%s' %('/mnt/titania-hecuba/mcintosh_lab/john/Data/notebooks', nb_name)
# symlink the two
import os
if not os.path.isdir(outdir_loc): os.makedirs(outdir_loc)
if not os.path.isdir(outdir): os.system('ln -s %s %s' %(outdir_loc,outdir))
# root path to tvb-library, tvb-data, etc. libraries.
tvb_folder = le['code_dir'] + '/libraries_of_mine/github'
# stuff for workdocs-cloudfiles
aws_key = 'drowssaperucesyreva'
aws_secret = '?teytidettopsuoyevah'
Importage
# Generic imports
import os,sys,glob,h5py,itertools,multiprocessing,\
numpy as np,pandas as pd
from datetime import datetime
from time import time
from scipy import optimize
from scipy.spatial import cKDTree
# Visualization stuff
%matplotlib inline
from matplotlib import pyplot as plt
from IPython.display import Image,display as d
import seaborn as sns
# Neuroimaging stuff
import nibabel as nib
from nilearn.plotting import plot_roi
from dipy.tracking.utils import apply_affine
# TVB stuff
nbso,nbse = sys.stdout,sys.stderr # hack part 1/2 to keep output printing properly
sys.path.append(tvb_folder +'/tvb-library')
sys.path.append(tvb_folder + '/tvb-data')
from tvb.simulator.lab import *
sys.stdout,sys.stderr = nbso,nbse # ...hack part 2/2
from tvb.simulator.plot.tools import plot_surface_mpl,plot_surface_mpl_mv
from copy import copy # seems that this needs to be done last otherwise it gets
# replaced with numpy copy function
# aws api and workdocs-cloudfiles stuff
sys.path.append(le['ipynb_workdocs_dir'])
from new_utils import nb_fig,cloudfiles_nb
Initialize aws api and workdocs-cloudfiles folder
cnb = cloudfiles_nb('aws', [aws_key,aws_secret])
cnb.initialize_folder(nb_name)
Calico document tools
%%javascript
IPython.load_extensions('calico-spell-check',
'calico-document-tools', 'calico-cell-tools');
Go to output folder
os.chdir(outdir)
Ok, let's get cracking
Load data¶
The surface data we have are left and right hemisphere vertices, faces, and region mappings. The region data we have are weights and tract lengths matrices.
vtx_noinfl_lh = np.loadtxt('coord_left_noninflated.txt')[:,1:]
vtx_infl_lh = np.loadtxt('coord_left_inflated.txt')[:,1:]
tri_lh = np.loadtxt('tiles_left.txt').astype(int)
vtx_noinfl_rh = np.loadtxt('coord_right_noninflated.txt')[:,1:]
vtx_infl_rh = np.loadtxt('coord_right_inflated.txt')[:,1:]
tri_rh = np.loadtxt('tiles_right.txt').astype(int)
rm_lh = np.loadtxt('fv91_vertexMap_left.txt')
rm_rh = np.loadtxt('fv91_vertexMap_right.txt')
weights = np.loadtxt('weights.txt')
lengths = np.loadtxt('tract_lengths.txt')
Combine hemispheres
tri_lhrh = np.concatenate([tri_lh,tri_rh+vtx_infl_lh.shape[0]])
vtx_noinfl_lhrh = np.concatenate([vtx_noinfl_lh,vtx_noinfl_rh])
vtx_infl_lhrh = np.concatenate([vtx_infl_lh,vtx_infl_rh])
hemi = np.concatenate([np.zeros_like(vtx_infl_lh[:,0]).astype(bool),
np.ones_like(vtx_infl_rh[:,0]).astype(bool)])
Initial plots
figf = outdir + '/pics/bihemi_noinfl_surf_struct.png'
plot_surface_mpl_mv(vtx=vtx_noinfl_lhrh,tri=tri_lhrh,hemi=hemi,
reorient='fs',
shade_kwargs={'edgecolors': 'k', 'linewidth': 0.1})
plt.savefig(figf, bbox_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. figf = outdir + '/pics/bihemi_infl_surf_struct.png'
plot_surface_mpl_mv(vtx=vtx_infl_lhrh,tri=tri_lhrh,hemi=hemi,
reorient='fs',
shade_kwargs={'edgecolors': 'k', 'linewidth': 0.1})
plt.savefig(figf, bbox_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. The lateral and medial plots there are ok, but clearly something wrong with the combined superior view in the middle...
Fix inflated surface coordinates¶
So it turns out that the inflated surface coordinates are overlapping in the x plane; i.e. the left and right hemispheres are not separated.
The easy way to see this (which is how I found the problem) is with scatter plots
figf = outdir + '/pics/lhrh_infl_surf_scatter.png'
fig, ax = plt.subplots(ncols=3, figsize=(12,3))
ax[0].scatter(vtx_infl_lh[:,0][idx], vtx_infl_lh[:,1][idx],c='b',alpha=0.5)
ax[1].scatter(vtx_infl_lh[:,0][idx], vtx_infl_lh[:,2][idx],c='b',alpha=0.5)
ax[2].scatter(vtx_infl_lh[:,1][idx], vtx_infl_lh[:,2][idx],c='b',alpha=0.5)
ax[0].scatter(vtx_infl_rh[:,0][idx], vtx_infl_rh[:,1][idx],c='r',alpha=0.5)
ax[1].scatter(vtx_infl_rh[:,0][idx], vtx_infl_rh[:,2][idx],c='r',alpha=0.5)
ax[2].scatter(vtx_infl_rh[:,1][idx], vtx_infl_rh[:,2][idx],c='r',alpha=0.5)
plt.savefig(figf, bbox_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. Compare this to the non-inflated surfaces
figf = outdir + '/pics/lhrh_noinfl_surf_scatter.png'
fig, ax = plt.subplots(ncols=3, figsize=(12,3))
ax[0].scatter(vtx_noinfl_lh[:,0][idx], vtx_noinfl_lh[:,1][idx],c='b',alpha=0.5)
ax[1].scatter(vtx_noinfl_lh[:,0][idx], vtx_noinfl_lh[:,2][idx],c='b',alpha=0.5)
ax[2].scatter(vtx_noinfl_lh[:,1][idx], vtx_noinfl_lh[:,2][idx],c='b',alpha=0.5)
ax[0].scatter(vtx_noinfl_rh[:,0][idx], vtx_noinfl_rh[:,1][idx],c='r',alpha=0.5)
ax[1].scatter(vtx_noinfl_rh[:,0][idx], vtx_noinfl_rh[:,2][idx],c='r',alpha=0.5)
ax[2].scatter(vtx_noinfl_rh[:,1][idx], vtx_noinfl_rh[:,2][idx],c='r',alpha=0.5)
plt.savefig(figf, bbox_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. We can fix this by shifting the x coordinates out to the left for the left hemisphere and out to the right for the right hemisphere:
vtx_infl_shift_lh = vtx_infl_lh.copy()
vtx_infl_shift_lh[:,0] -= vtx_infl_lh[:,0].max()
vtx_infl_shift_lh[:,0] -= 1.
vtx_infl_shift_rh = vtx_infl_rh.copy()
vtx_infl_shift_rh[:,0] += np.abs(vtx_infl_rh[:,0].min())
vtx_infl_shift_rh[:,0] += 1.
vtx_infl_shift_lhrh = np.concatenate([vtx_infl_shift_lh,vtx_infl_shift_rh])
figf = outdir + '/pics/lhrh_infl_shift_surf_scatter.png'
fig, ax = plt.subplots(ncols=3, figsize=(12,3))
ax[0].scatter(vtx_infl_shift_lh[:,0][idx], vtx_infl_shift_lh[:,1][idx],c='b',alpha=0.5)
ax[1].scatter(vtx_infl_shift_lh[:,0][idx], vtx_infl_shift_lh[:,2][idx],c='b',alpha=0.5)
ax[2].scatter(vtx_infl_shift_lh[:,1][idx], vtx_infl_shift_lh[:,2][idx],c='b',alpha=0.5)
ax[0].scatter(vtx_infl_shift_rh[:,0][idx], vtx_infl_shift_rh[:,1][idx],c='r',alpha=0.5)
ax[1].scatter(vtx_infl_shift_rh[:,0][idx], vtx_infl_shift_rh[:,2][idx],c='r',alpha=0.5)
ax[2].scatter(vtx_infl_shift_rh[:,1][idx], vtx_infl_shift_rh[:,2][idx],c='r',alpha=0.5)
plt.savefig(figf, bboc_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. And now the surface plots look fine as well
figf = outdir + '/pics/bihemi_infl_shift_surf_struct.png'
fig, ax = plt.subplots(ncols=3, figsize=(12,3))
plot_surface_mpl_mv(vtx=vtx_infl_shift_lhrh,tri=tri_lhrh,hemi=hemi,reorient='fs',
shade_kwargs={'edgecolors': 'k', 'linewidth': 0.1})
plt.savefig(figf, bboc_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. Plotting activity from a region mapping¶
Now we want to plot data from a set of regions that are linked to our surface vertices with a region mapping
The data we have from Gleb contain region mappings for the left and right hemispheres separately. We want to use these separately, and also together for combined left + right hemisphere plots.
A slight quick of these region mappings is that they contain a code (0) that indicates surface vertices that don't correspond to any of our brain regions. This makes things a little fiddly, as we'll see; but not too problematic.
We deal with this by appending a 'dummy' entry to then end of our data vectors, and changing the 0s in the region mapping codes to match this new dummy entry index.
We do this separately for each hemisphere and for the two hemispheres combined, so we end up with three 'augmented' data vectors.
(We're actually getting data from a connectivity weights/tract lengths in this case, so we have three augmented matrices. But we are only ever using one column at a time of these matrices.
nph = 77 # number of nodes per hemiphere
# Make augmented weights and tract lengths matrices
weights_aug_lh = np.zeros([nph+1,nph+1])
weights_aug_lh[:-1,:-1] = weights[:nph,:nph].copy()
weights_aug_rh = np.zeros([nph+1,nph+1])
weights_aug_rh[:-1,:-1] = weights[nph:,nph:].copy()
weights_aug_lhrh = np.zeros([nph*2+1,nph*2+1])
weights_aug_lhrh[:-1,:-1] = weights.copy()
lengths_aug_lh = np.zeros([nph+1,nph+1])
lengths_aug_lh[:-1,:-1] = weights[:nph,:nph].copy()
lengths_aug_rh = np.zeros([nph+1,nph+1])
lengths_aug_rh[:-1,:-1] = lengths[nph:,nph:].copy()
lengths_aug_lhrh = np.zeros([nph*2+1,nph*2+1])
lengths_aug_lhrh[:-1,:-1] = lengths.copy()
# Make augmented region mappings
rm_aug_lh = rm_lh.copy()
rm_aug_lh[rm_aug_lh==0] = 9999 # switch 0 code to this temporarily
rm_aug_lh-=1 # convert all indexes from matlab to python (0-based)
rm_aug_lh-=100 # indices for these should be 0:npm
rm_aug_lh = rm_aug_lh.astype(int)
rm_aug_rh = rm_rh.copy()
rm_aug_rh[rm_aug_rh==0] = 9999
rm_aug_rh-=1
rm_aug_rh = rm_aug_rh.astype(int)
rm_aug_lhrh = np.concatenate([rm_aug_lh,rm_aug_rh+77])
rm_aug_lh[rm_aug_lh>9000] = 78 # dummy entry index for lh data
rm_aug_rh[rm_aug_rh>9000] = 78 # dummy entry index for rh data
rm_aug_lhrh[rm_aug_lhrh>9000] = 154 # dummy entry index for bihemi data
So these are the augmented bihemispheric matrices.
figf = outdir + '/pics/weights_aug_bihemi.png'
fig, ax = plt.subplots(ncols=2,figsize=(12,3))
sns.heatmap(weights_aug_lhrh, xticklabels='',yticklabels='',ax=ax[0])
sns.heatmap(lengths_aug_lhrh, xticklabels='',yticklabels='',ax=ax[1])
plt.savefig(figf, bbox_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. Now we will do some surface plots for one columns of these matrices. The column is the second column (chosen at random), which corresponds to area 10 (frontal cortex).
roi_ind = 1
figf = outdir + '/pics/bihemi_infl_surf_weightsauglhrh_roi1.png'
dat=weights_aug_lhrh[:,roi_ind].copy()
# set the weight value for the roi to a high value so we can see the region
# prominently
dat[1] = 5
plot_surface_mpl_mv(vtx=vtx_infl_shift_lhrh,tri=tri_lhrh,hemi=hemi,reorient='fs',
data=dat,rm=rm_aug_lhrh)
plt.savefig(figf, bbox_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. So the red blob in the above figure is the actual region that the connectivity weights map corresponds to. Now plotting the distances for that same region to all other regions:
figf = outdir + '/pics/bihemi_infl_surf_lengthsauglhrh_roi1.png'
dat=lengths_aug_lhrh[:,roi_ind].copy()
# (won't bother changing the value for roi1 this time...)
plot_surface_mpl_mv(vtx=vtx_infl_shift_lhrh,tri=tri_lhrh,hemi=hemi,reorient='fs',
data=dat,rm=rm_aug_lhrh)
plt.savefig(figf, bbox_inches='tight')
plt.close()
# Upload to cloud
cap = ''
label = ""
im = nb_fig(figf,label,cap,cnb,
size=(1000,600),show_fignum=False)
d(im)
. Looking good to me :)