In [ ]:
try:
%load_ext autoreload
%autoreload 2
%reset
except:
print 'NOT IPYTHON'
from __future__ import division
import os
os.environ['MKL_NUM_THREADS'] = '1'
os.environ['NUMEXPR_NUM_THREADS'] = '1'
os.environ['OMP_NUM_THREADS'] = '1'
import sys
import glob
import numpy as np
import pandas as pd
import seaborn as sns
import scipy.stats as stats
import scipy.io as io
import h5py
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import rcParams
echobase_path = '/Users/akhambhati/Developer/hoth_research/Echobase'
#echobase_path = '/data/jag/akhambhati/hoth_research/Echobase'
sys.path.append(echobase_path)
import Echobase
convert_conn_vec_to_adj_matr = Echobase.Network.Transforms.configuration.convert_conn_vec_to_adj_matr
convert_adj_matr_to_cfg_matr = Echobase.Network.Transforms.configuration.convert_adj_matr_to_cfg_matr
subgraph = Echobase.Network.Partitioning.Subgraph
rcParams = Echobase.Plotting.fig_format.update_rcparams(rcParams)
path_Remotes = '/Users/akhambhati/Remotes'
#path_Remotes = '/data/jag/bassett-lab/akhambhati'
path_CoreData = path_Remotes + '/CORE.fMRI_cogcontrol.medaglia'
path_PeriphData = path_Remotes + '/RSRCH.NMF_CogControl'
path_InpData = path_PeriphData + '/e01-FuncNetw'
path_ExpData = path_PeriphData + '/e02-FuncSubg'
path_Figures = './e02-Figures'
for path in [path_CoreData, path_PeriphData, path_ExpData, path_Figures]:
if not os.path.exists(path):
print('Path: {}, does not exist'.format(path))
os.makedirs(path)
In [ ]:
# Load configuration matrix
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_InpData))
cfg_matr = df['cfg_matr']
cfg_obs_lut = df['cfg_obs_lut']
n_subj = len(df['cfg_key_label'][()]['Subject_ID'])
# Generate folds
n_fold = 4
n_obs_per_fold = np.size(cfg_obs_lut) / n_fold
assert n_obs_per_fold.is_integer()
n_obs_per_fold = int(n_obs_per_fold)
rand_obs = np.random.permutation(np.size(cfg_obs_lut))
fold_obs = rand_obs.reshape(n_fold, n_obs_per_fold)
fold_list = [list(ff) for ff in fold_obs]
# Cross-Validation Progress
str_path = '{}/NMF_CrossValidation.Progress.txt'.format(path_ExpData)
if os.path.exists(str_path):
os.remove(str_path)
# Get parameter search space
param_list = Echobase.Network.Partitioning.Subgraph.optimize_nmf.gen_random_sampling_paramset(
rank_range=(3, 51),
alpha_range=(1e-2, 5.0),
beta_range=(1e-2, 5.0),
n_param=1000,
fold_list=fold_list,
str_path=str_path)
# Save param_list for sge run
np.savez('{}/NMF_CrossValidation.Param_List.npz'.format(path_ExpData),
param_list=param_list)
In [ ]:
# Map NMF xval to all the parameter sets
job_str = './NMF_xval.py {} {} {}'.format(echobase_path, path_InpData, path_ExpData)
qsub_str = 'qsub -t 1-{} {}'.format(len(param_list), job_str)
os.chdir('./e02-SGE_Scripts/')
!sh {qsub_str}
os.chdir('../')
In [ ]:
# Reduce NMF xval output to a qmeas_list
path_xval_out = glob.glob('{}/NMF_CrossValidation.Param.*.npz'.format(path_ExpData))
qmeas_list = [np.load(pth)['qmeas_dict'][()] for pth in path_xval_out]
In [ ]:
param_list = np.load('{}/NMF_CrossValidation.Param_List.npz'.format(path_ExpData))['param_list']
all_param, opt_param = Echobase.Network.Partitioning.Subgraph.optimize_nmf.find_optimum_xval_paramset(param_list, qmeas_list, search_pct=5.0)
n_bin = 100
srch_pct = 25.0
for param in ['rank', 'alpha', 'beta']:
err_mean = stats.binned_statistic(all_param[param], all_param[qmeas],
statistic=np.nanmean, bins=n_bin)
err_std = stats.binned_statistic(all_param[param], all_param[qmeas],
statistic=np.nanstd, bins=n_bin)
mean_ix = np.flatnonzero(err_mean[0] < np.nanpercentile(err_mean[0], srch_pct))
std_ix = np.flatnonzero(err_std[0] < np.nanpercentile(err_std[0], srch_pct))
opt_param[param] = np.nanmean(err_mean[1][np.intersect1d(mean_ix, std_ix)])
if param == 'rank':
opt_param[param] = int(np.round(opt_param[param]))
print
print('Optimal Rank: {}'.format(opt_param['rank']))
print('Optimal Alpha: {}'.format(opt_param['alpha']))
print('Optimal Beta: {}'.format(opt_param['beta']))
np.savez('{}/NMF_CrossValidation.Optimal_Param.npz'.format(path_ExpData),
opt_param=opt_param,
all_param=all_param)
In [ ]:
### Amass all parameters and compute optima based on discussion in (Ng (1997). ICML)
opt_dict = np.load('{}/NMF_CrossValidation.Optimal_Param.npz'.format(path_ExpData))
opt_param = opt_dict['opt_param'][()]
all_param = opt_dict['all_param'][()]
# Generate quality measure plots
for qmeas in ['test_error']:
for param in ['rank', 'alpha', 'beta']:
ax_jp = sns.jointplot(all_param[param], all_param[qmeas], kind='kde',
space=0, n_levels=60, shade_lowest=False)
ax = ax_jp.ax_joint
ax.plot([opt_param[param], opt_param[param]],
[ax.get_ylim()[0], ax.get_ylim()[1]],
lw=1.0, alpha=0.75, linestyle='--')
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xlabel(param)
ax.set_ylabel(qmeas)
plt.savefig('{}/NMF_Optimization.{}.{}.svg'.format(path_Figures, param, qmeas))
plt.show()
plt.close()
In [ ]:
# Map NMF consensus to all the parameter sets
n_opt = 1000
job_str = './NMF_consensus_map.py {} {} {}'.format(echobase_path, path_InpData, path_ExpData)
qsub_str = 'qsub -t 1-{} {}'.format(n_opt, job_str)
os.chdir('./e02-SGE_Scripts/')
!sh {qsub_str}
os.chdir('../')
In [ ]:
# Reduce NMF consensus from all seed subgraphs
job_str = './NMF_consensus_reduce.py {} {} {}'.format(echobase_path, path_InpData, path_ExpData)
qsub_str = 'qsub {}'.format(job_str)
os.chdir('./e02-SGE_Scripts/')
!sh {qsub_str}
os.chdir('../')
In [ ]:
# Map NMF surrogate
n_opt = 1000
job_str = './NMF_surrogate_map.py {} {} {}'.format(echobase_path, path_InpData, path_ExpData)
qsub_str = 'qsub -t 1-{} {}'.format(n_opt, job_str)
os.chdir('./e02-SGE_Scripts/')
!sh {qsub_str}
os.chdir('../')
"""
for nn in xrange(n_opt):
if not os.path.exists('{}/NMF_Surrogate.Param.{}.npz'.format(path_ExpData, nn)):
qsub_str = 'qsub -t {}-{} {}'.format(nn+1, nn+1, job_str)
os.chdir('./e02-SGE_Scripts/')
!sh {qsub_str}
os.chdir('../')
"""
In [ ]:
%matplotlib inline
# Load the consensus data
data = np.load("{}/NMF_Consensus.Param.All.npz".format(path_ExpData), mmap_mode='r')
#data = np.load("{}/NMF_Surrogate.Param.50.npz".format(path_ExpData), mmap_mode='r')
fac_subnet = data['fac_subnet']
fac_coef = data['fac_coef']
n_fac = fac_subnet.shape[0]
n_conn = fac_subnet.shape[1]
n_win = fac_coef.shape[1]
# Plot subgraph matrix
plt.figure()
ax = plt.subplot(111)
mat = ax.matshow(fac_subnet.T, aspect=float(n_fac)/n_conn, cmap='rainbow')
#plt.colorbar(mat, ax=ax)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
#ax.set_xticks(np.linspace(0, 80, 5))
ax.set_ylabel('Functional Interactions')
ax.set_xlabel('Subgraphs')
plt.savefig('{}/Subgraph-Cfg_Matrix.svg'.format(path_Figures))
plt.show()
plt.close()
# Plot subgraph adjacency
plt.figure()
n_row = np.floor(np.sqrt(n_fac))
n_col = np.ceil(n_fac / n_row)
for ii, subg in enumerate(fac_subnet):
adj = convert_conn_vec_to_adj_matr(subg)
ax = plt.subplot(n_row, n_col, ii+1)
mat = ax.matshow(adj, cmap='rainbow') #, vmin=0, vmax=1)
#plt.colorbar(mat, ax=ax)
ax.set_axis_off()
plt.savefig('{}/Subgraph-Adj_Matrices.svg'.format(path_Figures))
plt.show()
plt.close()
# Plot Coefficients
plt.figure()
ax = plt.subplot(111)
fac_coef = fac_coef.T
norm_fac = fac_coef - fac_coef.mean(axis=0)
for ff in xrange(n_fac):
ax.plot(ff + norm_fac[:, ff] / (3*np.std(norm_fac[:, ff])), color=[66/256., 152/256., 221./256])
# Axis Settings
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_ylim([-1, n_fac+1])
ax.set_xlim([0, int(n_win/28*3)])
ax.set_ylabel('Subgraphs')
ax.set_xlabel('Time Windows')
plt.savefig('{}/Subgraph-Coefs.svg'.format(path_Figures))
plt.show()
plt.close()
In [ ]:
# Map NMF consensus to all the parameter sets
n_opt = 1000
job_str = './NMF_consensus_map-split_half.py {} {} {}'.format(echobase_path, path_InpData, path_ExpData)
qsub_str = 'qsub -t 1-{} {}'.format(n_opt, job_str)
os.chdir('./e02-SGE_Scripts/')
!sh {qsub_str}
os.chdir('../')
In [ ]:
# Reduce NMF consensus from all seed subgraphs
job_str = './NMF_consensus_reduce-split_half.py {} {} {}'.format(echobase_path, path_InpData, path_ExpData)
qsub_str = 'qsub {}'.format(job_str)
os.chdir('./e02-SGE_Scripts/')
!sh {qsub_str}
os.chdir('../')
In [ ]:
import scipy.optimize as sciopt
data_A = np.load('{}/NMF_Consensus.Param.All.npz'.format(path_ExpData),
mmap_mode='r')
fac_subnet_A = data_A['fac_subnet']
data_B = np.load('{}/NMF_Consensus.Param.A.All.npz'.format(path_ExpData),
mmap_mode='r')
fac_subnet_B = data_B['fac_subnet']
n_fac = fac_subnet_A.shape[0]
cost_matrix = np.zeros((n_fac, n_fac))
for fac_ii in xrange(fac_subnet_A.shape[0]):
for fac_jj in xrange(fac_subnet_B.shape[0]):
cost_matrix[fac_ii, fac_jj] = np.linalg.norm(fac_subnet_A[fac_ii] - fac_subnet_B[fac_jj])
old_A_ii, new_A_ii = sciopt.linear_sum_assignment(cost_matrix)
fac_subnet_A[new_A_ii, :] = fac_subnet_A[old_A_ii, :]
true_rho = np.array([stats.pearsonr(fac_subnet_A[fac_i, :], fac_subnet_B[fac_i, :])
for fac_i in xrange(n_fac)])
true_rho, true_pv = true_rho[:, 0], true_rho[:, 1]
null_rho = []
for fac_i in xrange(n_fac):
for fac_j in xrange(n_fac):
if fac_i == fac_j:
continue
null_rho.append(stats.pearsonr(fac_subnet_A[fac_i, :], fac_subnet_B[fac_j, :])[0])
null_rho = np.array(null_rho)
### Plot Test-Retest Reliability
plt.figure(dpi=300)
ax = plt.subplot(111)
clr = []
all_pv = []
for rho in np.sort(true_rho)[::-1]:
all_pv.append(np.mean(null_rho > rho))
for is_sig in Echobase.Statistics.FDR.fdr.bhp(all_pv):
if is_sig:
clr.append('b')
else:
clr.append('k')
ax.bar(xrange(len(true_rho)), np.sort(true_rho)[::-1], color=clr, lw=0)
ax.hlines(np.percentile(null_rho, 95), -1, n_fac, color='r', lw=1.0, linestyle='--');
ax.hlines(np.percentile(null_rho, 5), -1, n_fac, color='r', lw=1.0, linestyle='--');
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xlabel('Subgraph Pairs')
ax.set_ylabel('Pearson r Between Subgraphs')
ax.set_xlim([-1, n_fac])
ax.set_ylim([-0.5, 1.0])
plt.savefig('{}/Test_Retest.svg'.format(path_Figures))
plt.show()
In [ ]:
# Grab the subgraphs and expression from consensus NMF
df_nmf = np.load("{}/NMF_Consensus.Param.All.npz".format(path_ExpData),
mmap_mode='r')
fac_subnet = df_nmf['fac_subnet']
fac_coef = df_nmf['fac_coef']
n_fac = fac_subnet.shape[0]
n_conn = fac_subnet.shape[1]
n_node = np.int(np.ceil(np.sqrt(n_conn*2)))
n_obs = fac_coef.shape[1]
# Retrieve the configuration matrix
df_cfg = np.load('{}/Population.Configuration_Matrix.npz'.format(path_InpData))
cfg_obs_lut = np.array(df_cfg['cfg_obs_lut'], dtype=np.int)
# Retrieve the Yeo Systems Assignments for Lausanne 125
df_to_yeo = np.load('{}/Lausanne125_to_Yeo.npz'.format(path_InpData))
n_laus = len(df_to_yeo['yeo_lbl'])
n_yeo = len(df_to_yeo['yeo_name'])
In [ ]:
sys_subgraph_path = '{}/Yeo_Subgraph.All.npz'.format(path_ExpData)
if os.path.exists(sys_subgraph_path):
df_subg = np.load('{}/Yeo_Subgraph.All.npz'.format(path_ExpData))
else:
system_subgraph = []
for fac_ii in xrange(n_fac):
print('Processing: {} of {}'.format(fac_ii+1, n_fac))
sel_subnet = fac_subnet[fac_ii, :]
sel_coef = fac_coef[fac_ii, :]
adj_roi = convert_conn_vec_to_adj_matr(sel_subnet)
# Compute Brain System Adjacency Matrices
n_yeo = len(df_to_yeo['yeo_name'])
adj_yeo = np.zeros((n_yeo, n_yeo))
# Permutation
n_perm = 10000
alpha = 0.05 / len(df_to_yeo['yeo_triu'][0])
sel_subnet_null = np.array([np.random.permutation(sel_subnet)
for n_ii in xrange(n_perm)])
adj_yeo_null = np.zeros((n_perm, n_yeo, n_yeo))
# Mean inter/intra system edge wt
for sys_ix, sys_iy in zip(*df_to_yeo['yeo_triu']):
sys1 = df_to_yeo['yeo_name'][sys_ix]
sys2 = df_to_yeo['yeo_name'][sys_iy]
sys1_ix = np.flatnonzero(df_to_yeo['yeo_lbl'][df_to_yeo['laus_triu'][0]] == sys1)
sys1_iy = np.flatnonzero(df_to_yeo['yeo_lbl'][df_to_yeo['laus_triu'][1]] == sys1)
sys2_ix = np.flatnonzero(df_to_yeo['yeo_lbl'][df_to_yeo['laus_triu'][0]] == sys2)
sys2_iy = np.flatnonzero(df_to_yeo['yeo_lbl'][df_to_yeo['laus_triu'][1]] == sys2)
inter_sys_ii = np.unique(np.concatenate((np.intersect1d(sys1_ix, sys2_iy),
np.intersect1d(sys1_iy, sys2_ix))))
# Populate a full adjacency matrix
adj_yeo[sys_ix, sys_iy] = np.mean(sel_subnet[inter_sys_ii])
adj_yeo[sys_iy, sys_ix] = np.mean(sel_subnet[inter_sys_ii])
adj_yeo_null[:, sys_ix, sys_iy] = np.mean(sel_subnet_null[:, inter_sys_ii], axis=1)
adj_yeo_null[:, sys_iy, sys_ix] = np.mean(sel_subnet_null[:, inter_sys_ii], axis=1)
# Compute core-periphery scores
intra_sys = np.diag(adj_yeo)
inter_sys = np.sum(np.triu(adj_yeo, k=1) +
np.triu(adj_yeo, k=1).T, axis=1) / (n_yeo-1)
# Compute core-periphery null scores
null_intra_sys = np.array([np.diag(aa) for aa in adj_yeo_null])
null_inter_sys = np.array([np.sum(np.triu(aa, k=1) +
np.triu(aa, k=1).T, axis=1) / (n_yeo-1)
for aa in adj_yeo_null])
# Threshold
adj_yeo[(np.mean(adj_yeo_null > adj_yeo, axis=0) >= alpha)] = 0
# Generate subgraph dictionary
system_subgraph.append({'Subgraph_ID': fac_ii+1,
'subnet_yeo': adj_yeo,
'subnet_roi': adj_roi[df_to_yeo['sort_laus_to_yeo'], :][:, df_to_yeo['sort_laus_to_yeo']],
'intra_sys': intra_sys,
'inter_sys': inter_sys,
'null_intra_sys': null_intra_sys,
'null_inter_sys': null_inter_sys,
'expr_coef': sel_coef})
np.savez('{}/Yeo_Subgraph.All.npz'.format(path_ExpData),
system_subgraph=system_subgraph)
In [ ]:
# Plot the subgraphs
plt.figure(figsize=(5,5), dpi=300);
for ii, fac_ii in enumerate(xrange(n_fac)): #enumerate(sort_fac):
sel_fac_subnet = fac_subnet_B[fac_ii, :]
adj = convert_conn_vec_to_adj_matr(sel_fac_subnet)
adj_yeo = adj[df_to_yeo['sort_laus_to_yeo'], :][:, df_to_yeo['sort_laus_to_yeo']]
# Plot
ax = plt.subplot(4, 4, ii+1)
mat = ax.matshow(adj_yeo,
cmap='magma', vmin=1.1*sel_fac_subnet.min());
#plt.colorbar(mat, ax=ax, fraction=0.046, pad=0.04)
for xx in df_to_yeo['yeo_adj_demarc']:
ax.vlines(xx, 0, n_laus, color='w', lw=0.25)
ax.hlines(xx, 0, n_laus, color='w', lw=0.25)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_tick_params(width=0)
ax.xaxis.set_tick_params(width=0)
ax.grid(False)
ax.tick_params(axis='both', which='major', pad=-3)
ax.set_xticks((df_to_yeo['yeo_adj_demarc'][:-1] + (np.diff(df_to_yeo['yeo_adj_demarc']) * 0.5)));
ax.set_xticklabels(df_to_yeo['yeo_name'], fontsize=3.0, rotation=45)
ax.set_yticks((df_to_yeo['yeo_adj_demarc'][:-1] + (np.diff(df_to_yeo['yeo_adj_demarc']) * 0.5)));
ax.set_yticklabels(df_to_yeo['yeo_name'], fontsize=3.0, rotation=45)
plt.show()
In [ ]:
pixmap_path = '{}/brain_system_pixmap.npz'.format(path_Figures)
if os.path.exists(pixmap_path):
brain_system_pixmap = np.load(pixmap_path)
else:
from mayavi import mlab
import nibabel as nib
brain_system_pixmap = {}
brain_system_rgbcol = {}
sys_scalar = [15, 18, 6, 0, 3, 4, 2, 8, 10]
view_angle = {'Sag_PA': [0.0, 90.0],
'Sag_AP': [180.0, 90.0]}
# Get the pial surface recons
pial_hemi = {'LH': {},
'RH': {}}
pial_hemi['LH']['vert'], pial_hemi['LH']['tria'] = nib.freesurfer.io.read_geometry('{}/BrainRenderSubject15/surf/lh.pial'.format(path_CoreData))
pial_hemi['RH']['vert'], pial_hemi['RH']['tria'] = nib.freesurfer.io.read_geometry('{}/BrainRenderSubject15/surf/rh.pial'.format(path_CoreData))
# Get the Lausanne label files for each ROI
label_files = []
for roi in df_to_yeo['df_laus_yeo']:
laus_lbl = roi[1].lower()
hemi = roi[2].lower()
# Parse the atlas name and find the label file if it exists
lbl_file = '{}.{}.label'.format(hemi, laus_lbl)
lbl_file = lbl_file.replace(' ', '')
label_files.append('{}/BrainRenderSubject15/label/regenerated_{}_125/{}'.format(path_CoreData, hemi, lbl_file))
# Iterate over hemisphere of the pial surface
for hemi in pial_hemi.keys():
n_vert = len(pial_hemi[hemi]['vert'])
# Iterate over brain system
for sys_id, sys_lbl in enumerate(df_to_yeo['yeo_name']):
print(sys_lbl)
sys_ix = np.flatnonzero(df_to_yeo['yeo_lbl'] == sys_lbl)
# Find the label file for each ROI and get vertices
if sys_lbl == 'SUB':
pial_scalars = sys_scalar[sys_id]*np.ones(n_vert)
else:
pial_scalars = 15*np.ones(n_vert)
for roi_ix, (roi, lbl_file) in enumerate(zip(df_to_yeo['df_laus_yeo'], label_files)):
if roi[2] != hemi:
continue
if not os.path.exists(lbl_file):
continue
# Load the file and add scalar to the vertices
parc_lbl = nib.freesurfer.io.read_label(lbl_file)
if roi_ix in sys_ix:
pial_scalars[parc_lbl] = sys_scalar[sys_id]
else:
pial_scalars[parc_lbl] = 15
# Plot the colored Brain System
fig = mlab.figure(bgcolor=(1.0, 1.0, 1.0))
src = mlab.pipeline.triangular_mesh_source(pial_hemi[hemi]['vert'][:,0],
pial_hemi[hemi]['vert'][:,1],
pial_hemi[hemi]['vert'][:,2],
pial_hemi[hemi]['tria'], scalars=pial_scalars, opacity=0.75, figure=fig)
norms = mlab.pipeline.poly_data_normals(src, figure=fig)
norms.filter.splitting = False
surf = mlab.pipeline.surface(norms, figure=fig)
surf.parent.scalar_lut_manager.set(lut_mode='Vega20', data_range=[0, 20], use_default_range=False)
lut = surf.module_manager.scalar_lut_manager.lut.table.to_array()
lut[188:213, 3] = 220
surf.module_manager.scalar_lut_manager.lut.table = lut
# Rotate the view and save a screenshot
pixmap = {}
for ang in view_angle.keys():
mlab.view(azimuth=view_angle[ang][0],
elevation=view_angle[ang][1])
pixmap[ang] = mlab.screenshot(mode='rgba')
mlab.close(all=True)
# Save to system pixmap dictionary
if not sys_lbl in brain_system_pixmap.keys():
brain_system_pixmap[sys_lbl] = {}
if not hemi in brain_system_pixmap[sys_lbl].keys():
brain_system_pixmap[sys_lbl][hemi] = pixmap
# Save RGB
brain_system_rgbcol[sys_lbl] = lut[(sys_scalar[sys_id]*lut.shape[0]/20.), :]
np.savez(pixmap_path,
brain_system_pixmap=brain_system_pixmap,
brain_system_rgbcol=brain_system_rgbcol)
In [ ]:
sys_pixmap = brain_system_pixmap['brain_system_pixmap'][()]
for sys_id in sys_pixmap.keys():
for hemi_id in ['RH']: #sys_pixmap[sys_id].keys():
for plane_id in sys_pixmap[sys_id][hemi_id].keys():
plt.figure(dpi=300.0)
ax = plt.subplot(111)
ax.imshow(sys_pixmap[sys_id][hemi_id][plane_id]);
ax.set_axis_off();
plt.savefig('{}/System_Pixmap.{}.{}.{}.svg'.format(path_Figures, sys_id, hemi_id, plane_id))
plt.close()
In [ ]:
from matplotlib.offsetbox import AnnotationBbox, OffsetImage, TextArea
from scipy import ndimage
# Assign color for each circle node based on associated brain system color
node_clr = np.array([brain_system_pixmap['brain_system_rgbcol'][()][sys_lbl][:3] / 255
for sys_lbl in df_to_yeo['yeo_lbl'][df_to_yeo['sort_laus_to_yeo']]])
# Arrange nodes around the circle
system_pos = []
node_rads = np.linspace(0, 2*np.pi - (2*np.pi/n_laus), n_laus)
for sys_ii, sys_nm in enumerate(df_to_yeo['yeo_name']):
sys_rad = np.mean(node_rads[df_to_yeo['yeo_lbl'] == sys_nm])
dd = 16
system_pos.append((sys_rad, dd))
# Render the circle plot
for fac_ii in xrange(n_fac):
sel_subnet = df_subg['system_subgraph'][fac_ii]['subnet_roi']
sys_con = convert_adj_matr_to_cfg_matr(np.expand_dims(sel_subnet, axis=0)).squeeze()
fig, ax = Echobase.Plotting.render_circle_connectivity.draw(conn_list=sys_con,
conn_pct=[99, 100],
conn_cmap='YlGnBu',
conn_linewidth=0.5,
node_color=node_clr)
fig.savefig('{}/Circle_Subgraph.{}.svg'.format(path_Figures, fac_ii+1))
plt.close()
In [ ]: