<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Assess network impact of lesion
Step1:  We now use the synthetic lesion constructed in the previous example in a ConWhAt lesion analysis.
Step2:  Take another quick look at this mask
Step3:  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.
Step4:  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
Step5:  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.
Step6:  Now, compute lesion overlap statistics.
Step7:  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
Step8:  Typically we will be mainly interested in two of these metric scores
Step9:  We can obtain these numbers as a 'modification matrix' (connectivity matrix)
Step10:  These two maps are, unsurprisingly, very similar
Step11:  (...with an alternative color scheme...)
Step12:  We can list directly the most affected (greatest % overlap) connections,
Step13:  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
Step14:  Now plotting on a glass brain
<ASSISTANT_TASK:>
Python Code:
# 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
Explanation: Assess network impact of lesion
End of explanation
lesion_file = 'synthetic_lesion_20mm_sphere_-46_-60_6.nii.gz' # we created this file from scratch in the previous example
Explanation: We now use the synthetic lesion constructed in the previous example in a ConWhAt lesion analysis.
End of explanation
lesion_img = nib.load(lesion_file)
plot_roi(lesion_file);
Explanation: Take another quick look at this mask:
End of explanation
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)
Explanation: 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.
End of explanation
cw_vca = VolConnAtlas(atlas_dir=atlas_dir)
Explanation: 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
End of explanation
idxs = 'all' # alternatively, something like: range(1,100), indicates the first 100 cnxns (rows in .vmfs)
Explanation: 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.
End of explanation
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)
Explanation: Now, compute lesion overlap statistics.
End of explanation
lo_df.head()
Explanation: 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:
End of explanation
lo_df[['TPR', 'corr_thrbin']].iloc[:10].T
Explanation: 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)
End of explanation
tpr_adj = nx.to_pandas_adjacency(lo_nx,weight='TPR')
cpr_adj = nx.to_pandas_adjacency(lo_nx,weight='corr_thrbin')
Explanation: We can obtain these numbers as a 'modification matrix' (connectivity matrix)
End of explanation
np.corrcoef(tpr_adj.values.ravel(), cpr_adj.values.ravel())
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]);
Explanation: These two maps are, unsurprisingly, very similar:
End of explanation
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]);
Explanation: (...with an alternative color scheme...)
End of explanation
cw_vca.vfms.loc[lo_df.index].head()
Explanation: We can list directly the most affected (greatest % overlap) connections,
End of explanation
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())
Explanation: 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:
End of explanation
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)
fig, ax = plt.subplots(figsize=(16,6))
plot_connectome(cpr_adj.values,ccs_arr,axes=ax)
Explanation: Now plotting on a glass brain:
End of explanation
<END_TASK>