In [1]:
# ConWhAt stuff
from conwhat import VolConnAtlas,StreamConnAtlas,VolTractAtlas,StreamTractAtlas
from conwhat.viz.volume import plot_vol_scatter
# Neuroimaging stuff
import nibabel as nib
from nilearn.plotting import (plot_stat_map,plot_surf_roi,plot_roi,
plot_connectome,find_xyz_cut_coords)
from nilearn.image import resample_to_img
# Viz stuff
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
# Generic stuff
import glob, numpy as np, pandas as pd, networkx as nx
from datetime import datetime
We now use the synthetic lesion constructed in the previous example in a ConWhAt lesion analysis.
In [2]:
lesion_file = 'synthetic_lesion_20mm_sphere_-46_-60_6.nii.gz' # we created this file from scratch in the previous example
Take another quick look at this mask:
In [3]:
lesion_img = nib.load(lesion_file)
plot_roi(lesion_file);
Since our lesion mask does not (by construction) have a huge amount of spatial detail, it makes sense to use one of the lower-resolution atlas. As one might expect, computation time is considerably faster for lower-resolution atlases.
In [4]:
cw_atlases_dir = '/global/scratch/hpc3230/Data/conwhat_atlases' # change this accordingly
atlas_name = 'CWL2k8Sc33Vol3d100s_v01'
atlas_dir = '%s/%s' %(cw_atlases_dir, atlas_name)
See the previous tutorial on 'exploring the conwhat atlases' for more info on how to examine the components of a given atlas in ConWhAt.
Initialize the atlas
In [5]:
cw_vca = VolConnAtlas(atlas_dir=atlas_dir)
Choose which connections to evaluate.
This is normally an array of numbers indexing entries in cw_vca.vfms.
Pre-defining connection subsets is a useful way of speeding up large analyses, especially if one is only interested in connections between specific sets of regions.
As we are using a relatively small atlas, and our lesion is not too extensive, we can assess all connections.
In [6]:
idxs = 'all' # alternatively, something like: range(1,100), indicates the first 100 cnxns (rows in .vmfs)
Now, compute lesion overlap statistics.
In [18]:
jlc_dir = '/global/scratch/hpc3230/joblib_cache_dir' # this is the cache dir where joblib writes temporary files
lo_df,lo_nx = cw_vca.compute_hit_stats(lesion_file,idxs,n_jobs=4,joblib_cache_dir=jlc_dir)
This takes about 20 minutes to run.
vca.compute_hit_stats() returns a pandas dataframe, lo_df, and a networkx object, lo_nx.
Both contain mostly the same information, which is sometimes more useful in one of these formats and sometimes in the other.
lo_df is a table, with rows corresponding to each connection, and columns for each of a wide set of statistical metrics for evaluating sensitivity and specificity of binary hit/miss data:
In [28]:
lo_df.head()
Out[28]:
Typically we will be mainly interested in two of these metric scores:
TPR - True positive (i.e. hit) rate: number of true positives, divided by number of true positives + number of false negatives
corr_thrbin - Pearson correlation between the lesion amge and the thresholded, binarized connectome edge image (group-level visitation map)
In [27]:
lo_df[['TPR', 'corr_thrbin']].iloc[:10].T
Out[27]:
We can obtain these numbers as a 'modification matrix' (connectivity matrix)
In [33]:
tpr_adj = nx.to_pandas_adjacency(lo_nx,weight='TPR')
cpr_adj = nx.to_pandas_adjacency(lo_nx,weight='corr_thrbin')
These two maps are, unsurprisingly, very similar:
In [104]:
np.corrcoef(tpr_adj.values.ravel(), cpr_adj.values.ravel())
Out[104]:
In [79]:
fig, ax = plt.subplots(ncols=2, figsize=(12,4))
sns.heatmap(tpr_adj,xticklabels='',yticklabels='',vmin=0,vmax=0.5,ax=ax[0]);
sns.heatmap(cpr_adj,xticklabels='',yticklabels='',vmin=0,vmax=0.5,ax=ax[1]);
(...with an alternative color scheme...)
In [70]:
fig, ax = plt.subplots(ncols=2, figsize=(12,4))
sns.heatmap(tpr_adj, xticklabels='',yticklabels='',cmap='Reds',
mask=tpr_adj.values==0,vmin=0,vmax=0.5,ax=ax[0]);
sns.heatmap(cpr_adj,xticklabels='',yticklabels='',cmap='Reds',
mask=cpr_adj.values==0,vmin=0,vmax=0.5,ax=ax[1]);
We can list directly the most affected (greatest % overlap) connections,
In [85]:
cw_vca.vfms.loc[lo_df.index].head()
Out[85]:
To plot the modification matrix information on a brain, we first need to some spatial locations to plot as nodes. For these, we calculate (an approprixation to) each atlas region's centriod location:
In [101]:
parc_img = cw_vca.region_nii
parc_dat = parc_img.get_data()
parc_vals = np.unique(parc_dat)[1:]
ccs = {roival: find_xyz_cut_coords(nib.Nifti1Image((dat==roival).astype(int),img.affine),
activation_threshold=0) for roival in roivals}
ccs_arr = np.array(ccs.values())
Now plotting on a glass brain:
In [ ]:
fig, ax = plt.subplots(figsize=(16,6))
plot_connectome(tpr_adj.values,ccs_arr,axes=ax,edge_threshold=0.2,colorbar=True,
edge_cmap='Reds',edge_vmin=0,edge_vmax=1.,
node_color='lightgrey',node_kwargs={'alpha': 0.4});
#edge_vmin=0,edge_vmax=1)
In [118]:
fig, ax = plt.subplots(figsize=(16,6))
plot_connectome(cpr_adj.values,ccs_arr,axes=ax)
Out[118]: