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
from statsmodels.sandbox.stats import multicomp
import scipy.stats as stats
import scipy.io as io
import h5py
import matplotlib
import matplotlib.pyplot as plt
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
rcParams = Echobase.Plotting.fig_format.update_rcparams(matplotlib.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_Netw = path_PeriphData + '/e01-FuncNetw'
path_InpData_Subg = path_PeriphData + '/e02-FuncSubg'
path_ExpData = path_PeriphData + '/e03-FuncSubg_Dynamics'
path_Figures = './e03-Figures/'
for path in [path_CoreData, path_PeriphData, path_ExpData]:
if not os.path.exists(path):
print('Path: {}, does not exist'.format(path))
os.makedirs(path)
In [ ]:
%matplotlib inline
# Load Functional Data
df_cfg = np.load('{}/Population.Configuration_Matrix.npz'.format(path_InpData_Netw))
df_subgraph = np.load('{}/Yeo_Subgraph.All.npz'.format(path_InpData_Subg))
df_to_yeo = np.load('{}/Lausanne125_to_Yeo.npz'.format(path_InpData_Netw))
df_nmf = np.load("{}/NMF_Consensus.Param.All.npz".format(path_InpData_Subg),
mmap_mode='r')
n_subj, _, _, _, _, n_block = df_cfg['cfg_obs_lut'].shape
n_fac = len(df_subgraph['system_subgraph'])
surr_coef = np.array([np.load(pth, mmap_mode='r')['fac_coef'][...]
for pth in glob.glob('{}/NMF_Surrogate.Param.*.npz'.format(path_InpData_Subg))])
surr_coef = surr_coef.reshape(-1, surr_coef.shape[-1])
# Load Behavioral Data
df_blk = io.loadmat('{}/BlockwiseDataCorrectTrialsOnly.mat'.format(path_CoreData))
bad_subj_ix = [1, 6]
good_subj_ix = np.setdiff1d(np.arange(n_subj+2), bad_subj_ix)
"""
'high control accuracy', 'low control accuracy', 'high control mean RT',
'high control median RT', 'low control mean RT', 'low control median RT'
"""
df_perf = {'Stroop': {'lo': {'accuracy': df_blk['StroopData'][good_subj_ix, 1, :],
'meanRT': df_blk['StroopData'][good_subj_ix, 4, :],
'medianRT': df_blk['StroopData'][good_subj_ix, 5, :]},
'hi': {'accuracy': df_blk['StroopData'][good_subj_ix, 0, :],
'meanRT': df_blk['StroopData'][good_subj_ix, 2, :],
'medianRT': df_blk['StroopData'][good_subj_ix, 3, :]}
},
'Navon' : {'lo': {'accuracy': df_blk['NavonData'][good_subj_ix, 1, :],
'meanRT': df_blk['NavonData'][good_subj_ix, 4, :],
'medianRT': df_blk['NavonData'][good_subj_ix, 5, :]},
'hi': {'accuracy': df_blk['NavonData'][good_subj_ix, 0, :],
'meanRT': df_blk['NavonData'][good_subj_ix, 2, :],
'medianRT': df_blk['NavonData'][good_subj_ix, 3, :]}
}
}
# Load Motion Data
df_motion = {'Stroop': io.loadmat('{}/StroopMove.mat'.format(path_CoreData))['move'][:, 0],
'Navon': io.loadmat('{}/NavonMove.mat'.format(path_CoreData))['move'][:, 0]}
In [ ]:
cfg_matr = df_cfg['cfg_matr'].copy()
#cfg_matr = (df_cfg['cfg_matr'].T * df_cfg['cfg_L2']).T
rank = 13
alpha = 1e-7
beta = 1e-7
n_fac = rank
# Grab the task ID of the current job (and the associated parameter dictionary)
fac_subnet = np.random.uniform(low=0, high=1.0, size=(rank, cfg_matr.shape[1]))
fac_coef = np.random.uniform(low=0, high=1.0, size=(rank, cfg_matr.shape[0]))
# Run NMF Algorithm
fac_subnet, fac_coef, err = Echobase.Network.Partitioning.Subgraph.nmf.snmf_bcd(
cfg_matr, alpha=alpha, beta=beta, fac_subnet_init=fac_subnet, fac_coef_init=fac_coef,
max_iter=20, sparse_dim='conn', verbose=True)
In [ ]:
import string
abcd = list(string.ascii_uppercase)
coef_ix = np.array(df_cfg['cfg_obs_lut'], dtype=int)
# Re-rank subgraphs based on positive/negative expression
del_expr_mean = []
del_expr_stdv = []
for fac_ii in xrange(n_fac):
sel_fac_coef = df_subgraph['system_subgraph'][fac_ii]['expr_coef'][coef_ix]
pos_expr = sel_fac_coef[:, :, :, :, 0, :]
neg_expr = sel_fac_coef[:, :, :, :, 1, :]
del_expr = (pos_expr-neg_expr).mean(axis=-1).mean(axis=-1).mean(axis=-1).mean(axis=-1)
del_expr_mean.append(del_expr.mean())
del_expr_stdv.append(del_expr.std() / np.sqrt(n_subj))
del_expr_mean = np.array(del_expr_mean)
del_expr_stdv = np.array(del_expr_stdv)
sort_fac = np.argsort(del_expr_mean)[::-1]
# Create a sorted dictionary
sort_fac_dict = {}
for ltr, fac_ii in zip(abcd, sort_fac):
sort_fac_dict[ltr] = fac_ii
# Plot distribution of mean relative expression
plt.figure(dpi=300)
ax = plt.subplot(111)
ax.bar(xrange(n_fac), del_expr_mean[sort_fac], yerr=del_expr_stdv[sort_fac], lw=0)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticks(np.arange(n_fac)+0.4);
ax.set_xticklabels(np.sort(sort_fac_dict.keys()));
ax.set_xlim([0, n_fac])
ax.set_xlabel('Subgraphs')
ax.set_ylabel('Mean Relative Expression')
plt.savefig('{}/Subgraph_Relative_Expression.svg'.format(path_Figures))
plt.show()
In [ ]:
fac_key = 'K'
fac_ii = sort_fac_dict[fac_key]
sel_fac_coef = df_subgraph['system_subgraph'][fac_ii]['expr_coef'][coef_ix]
pos_expr = sel_fac_coef[:, :, :, :, 0, :]
neg_expr = sel_fac_coef[:, :, :, :, 1, :]
# Plot distribution of mean relative expression
plt.figure(dpi=300)
ax = plt.subplot(111)
ax.plot(pos_expr[27, ...].reshape(-1), color='r')
ax.plot(neg_expr[27, ...].reshape(-1), color='b')
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xlabel('Blocks')
ax.set_ylabel('Expression')
plt.savefig('{}/Example_Pos_Neg.Subgraph_{}.svg'.format(path_Figures, fac_key))
plt.show()
plt.close()
In [ ]:
alpha = 0.05
motion_rv = []
motion_pv = []
motion_color = []
motion_fac = []
coef_ix = np.array(df_cfg['cfg_obs_lut'], dtype=int)
for fac_key in np.sort(sort_fac_dict.keys()):
fac_ii = sort_fac_dict[fac_key]
sel_fac_coef = df_subgraph['system_subgraph'][fac_ii]['expr_coef'][coef_ix]
fac_expr_subj = sel_fac_coef.reshape(n_subj, -1).mean(axis=-1)
rv, pv = stats.spearmanr(0.5*(df_motion['Stroop'] + df_motion['Navon']),
fac_expr_subj)
motion_rv.append(rv)
motion_pv.append(pv)
is_sig = multicomp.multipletests(motion_pv, alpha=alpha, method='fdr_bh')[0]
for sig_bool, fac_key in zip(is_sig, np.sort(sort_fac_dict.keys())):
if sig_bool:
motion_color.append('r')
motion_fac.append(fac_key)
else:
motion_color.append('k')
motion_rv = np.array(motion_rv)
motion_color = np.array(motion_color)
motion_fac = np.array(motion_fac)
# Identify Subgraphs that correlate with motion
plt.figure(dpi=300)
ax = plt.subplot(111)
ax.bar(xrange(n_fac), motion_rv,
lw=0, color=motion_color)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticks(np.arange(n_fac)+0.4);
ax.set_xticklabels(np.sort(sort_fac_dict.keys()));
ax.set_ylim([-1, 1])
ax.set_xlim([0, n_fac])
ax.set_xlabel('Subgraphs')
ax.set_ylabel('Spearman rho(Expression, Motion)')
plt.savefig('{}/Subgraph_Motion_Score.svg'.format(path_Figures))
plt.show()
In [ ]:
%matplotlib inline
n_laus = len(df_to_yeo['yeo_lbl'])
n_yeo = len(df_to_yeo['yeo_name'])
fsize = 5.5
fig = plt.figure(figsize=(6, 6), dpi=300)
for ii, fac_key in enumerate(np.sort(sort_fac_dict.keys())):
fac_ii = sort_fac_dict[fac_key]
sel_fac_subnet = df_subgraph['system_subgraph'][fac_ii]['subnet_roi']
fac_vec = convert_adj_matr_to_cfg_matr(np.expand_dims(sel_fac_subnet, axis=0))[0, :]
vmin, vmax = fac_vec.min(), fac_vec.max()
ax = fig.add_subplot(3, 4, ii+1)
mat = ax.matshow(sel_fac_subnet, cmap='magma', vmin=vmin, vmax=vmax)
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_title(fac_key)
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=5.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=5.0, rotation=45)
fig.tight_layout(pad=0.01, h_pad=0.01, w_pad=0.01)
fig.savefig('{}/Ranked_System_Subgraph.svg'.format(path_Figures))
fig.show()
In [ ]:
%matplotlib inline
coef_ix = np.array(df_cfg['cfg_obs_lut'], dtype=int)
mean_rel_expr = {}
corep_index = {}
n_laus = len(df_to_yeo['yeo_lbl'])
n_yeo = len(df_to_yeo['yeo_name'])
alpha = 0.05 / (2*n_yeo)
fsize = 5.5
fig = plt.figure(figsize=(6, 6), dpi=300)
for ii, fac_key in enumerate(np.sort(sort_fac_dict.keys())):
fac_ii = sort_fac_dict[fac_key]
sel_fac_coef = df_subgraph['system_subgraph'][fac_ii]['expr_coef'][coef_ix]
pos_expr = sel_fac_coef[:, :, :, :, 0, :]
neg_expr = sel_fac_coef[:, :, :, :, 1, :]
del_expr = (pos_expr-neg_expr).mean(axis=-1).mean(axis=-1).mean(axis=-1).mean(axis=-1).mean(axis=-1)
mean_rel_expr[fac_key] = del_expr
intra_sys = df_subgraph['system_subgraph'][fac_ii]['intra_sys']
inter_sys = df_subgraph['system_subgraph'][fac_ii]['inter_sys']
null_intra_sys = df_subgraph['system_subgraph'][fac_ii]['null_intra_sys']
null_inter_sys = df_subgraph['system_subgraph'][fac_ii]['null_inter_sys']
sig_intra_sys = np.mean(null_intra_sys > intra_sys, axis=0) < alpha
sig_inter_sys = np.mean(null_inter_sys > inter_sys, axis=0) < alpha
intra_sys_thresh = intra_sys.copy()
inter_sys_thresh = inter_sys.copy()
intra_sys_thresh[~sig_intra_sys] = 0
inter_sys_thresh[~sig_inter_sys] = 0
corep_index[fac_key] = np.nanmean((intra_sys_thresh - inter_sys_thresh) / \
(intra_sys_thresh + inter_sys_thresh))
intra_colors = np.array(['k' for xx in xrange(n_yeo)])
inter_colors = np.array(['k' for xx in xrange(n_yeo)])
intra_colors[sig_intra_sys] = 'r'
inter_colors[sig_inter_sys] = 'r'
ax = fig.add_subplot(3, 4, ii+1)
ax.barh(np.arange(n_yeo), intra_sys, color=intra_colors, lw=0)
ax.barh(np.arange(n_yeo), -1*inter_sys, color=inter_colors, lw=0)
ax.plot([0.0, 0.0], [0.0, n_yeo], color='r', lw=1.0)
for x_ii, (low_x, high_x) in enumerate(zip(np.percentile(null_intra_sys, 100-100*alpha, axis=0),
-1*np.percentile(null_inter_sys, 100-100*alpha, axis=0))):
ax.fill_between([low_x, high_x],
x_ii-0.1, x_ii+0.9,
color=[0.1, 0.1, 0.1], lw=0, alpha=0.25)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.tick_params(axis='both', which='major', pad=0)
ax.set_ylim([0, n_yeo])
ax.set_yticks(np.arange(n_yeo) + 0.4);
ax.set_yticklabels(df_to_yeo['yeo_name'], fontsize=5.0, rotation=45)
ax.set_xlim([-0.015, 0.015])
ax.set_xticks(np.linspace(-0.015, 0.015, 3));
ax.set_title(fac_key)
ax.invert_yaxis()
plt.savefig('{}/Ranked_System_Subgraph.CorePeriphery.svg'.format(path_Figures))
plt.show()
plt.figure(figsize=(2,2), dpi=300.0)
ax = plt.subplot(111)
print(stats.spearmanr(corep_index.values(), mean_rel_expr.values()))
mm, yy, _, _, _ = stats.linregress(corep_index.values(), mean_rel_expr.values())
ax.scatter(corep_index.values(), mean_rel_expr.values(),
lw=0, alpha=0.75, color=[0.25, 0.25, 0.25])
ax.plot(np.array([-1.0, 1.0]), mm*np.array([-1.0, 1.0])+yy, color='k')
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.tick_params(axis='both', which='major', pad=1)
ax.set_xlim([-1.1, 1.1])
ax.set_xticks(np.linspace(-1.0, 1.0, 5));
ax.set_xlabel('Core-Periphery')
min_val, max_val = -3.5, 8
ax.set_ylim([min_val, max_val])
ax.set_ylabel('Mean Relative Expression')
plt.savefig('{}/Relative_Expression.CorePeriphery.svg'.format(path_Figures))
plt.show()
In [ ]:
%matplotlib inline
coef_ix = np.array(df_cfg['cfg_obs_lut'], dtype=int)
alpha = 0.05
del_expr_mean_stroop = []
del_expr_stdv_stroop = []
del_expr_mean_navon = []
del_expr_stdv_navon = []
tv = []
pv = []
clr = []
for ii, fac_key in enumerate(np.sort(sort_fac_dict.keys())):
fac_ii = sort_fac_dict[fac_key]
sel_fac_coef = df_subgraph['system_subgraph'][fac_ii]['expr_coef'][coef_ix]
pos_expr = sel_fac_coef[:, :, :, :, 0, :]
neg_expr = sel_fac_coef[:, :, :, :, 1, :]
del_expr = (pos_expr-neg_expr).mean(axis=-1)[:, :, :, 1].mean(axis=-1)
del_expr_mean = del_expr.mean(axis=0)
del_expr_stdv = del_expr.std(axis=0) / np.sqrt(n_subj)
del_expr_mean_stroop.append(del_expr_mean[0])
del_expr_stdv_stroop.append(del_expr_stdv[0])
del_expr_mean_navon.append(del_expr_mean[1])
del_expr_stdv_navon.append(del_expr_stdv[1])
# Within-experiment Stats
tv_expr, pv_expr = stats.ttest_rel(*del_expr.T)
tv.append(tv_expr)
pv.append(pv_expr)
#
is_sig = multicomp.multipletests(pv, alpha=alpha, method='fdr_bh')[0]
for sig_bool in is_sig:
if sig_bool:
clr.append('r')
else:
clr.append([0.2, 0.2, 0.2])
print(stats.spearmanr(del_expr_mean_stroop, del_expr_mean_navon))
plt.figure(figsize=(3,3), dpi=300)
ax = plt.subplot(111)
for ii in xrange(n_fac):
xx = del_expr_mean_stroop[ii]
yy = del_expr_mean_navon[ii]
ax.scatter(xx, yy, lw=0, color=clr[ii], s=10.0)
ax.plot([xx, 1/2*(xx+yy)],
[yy, 1/2*(xx+yy)], color=clr[ii])
min_val, max_val = -3.5, 8
ax.plot([min_val, max_val],
[min_val, max_val], 'k', alpha=0.1)
ax.vlines(0, min_val, max_val, 'r')
ax.hlines(0, min_val, max_val, 'r')
ax.set_xlim([min_val, max_val])
ax.set_ylim([min_val, max_val])
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xlabel('Stroop: Relative Expression')
ax.set_ylabel('Navon: Relative Expression')
plt.savefig('{}/Stroop_vs_Navon.svg'.format(path_Figures))
plt.show()
plt.close()
In [ ]:
%matplotlib inline
coef_ix = np.array(df_cfg['cfg_obs_lut'], dtype=int)
for expr_ii, expr_id in enumerate(df_cfg['cfg_key_label'][()]['Experiment_ID']):
del_expr_mean_lo = []
del_expr_stdv_lo = []
del_expr_mean_hi = []
del_expr_stdv_hi = []
pv = []
tv = []
clr = []
for ii, fac_key in enumerate(np.sort(sort_fac_dict.keys())):
fac_ii = sort_fac_dict[fac_key]
sel_fac_coef = df_subgraph['system_subgraph'][fac_ii]['expr_coef'][coef_ix]
pos_expr = sel_fac_coef[:, expr_ii, :, :, :, :][:, :, :, 0, :][:, :, 1, :]
neg_expr = sel_fac_coef[:, expr_ii, :, :, :, :][:, :, :, 1, :][:, :, 1, :]
del_expr = (pos_expr-neg_expr).mean(axis=-1)
del_expr_mean = del_expr.mean(axis=0)
del_expr_stdv = del_expr.std(axis=0) / np.sqrt(n_subj)
del_expr_mean_lo.append(del_expr_mean[0])
del_expr_stdv_lo.append(del_expr_stdv[0])
del_expr_mean_hi.append(del_expr_mean[1])
del_expr_stdv_hi.append(del_expr_stdv[1])
# Within-experiment Stats
tv_expr, pv_expr = stats.ttest_rel(*del_expr.T)
pv.append(pv_expr)
print(tv_expr, pv_expr)
for is_sig in Echobase.Statistics.FDR.fdr.bhp(pv, alpha=0.05, dependent=True):
if is_sig:
clr.append('r')
else:
clr.append([0.2, 0.2, 0.2])
print(stats.spearmanr(del_expr_mean_lo, del_expr_mean_hi))
plt.figure(figsize=(3,3), dpi=300)
ax = plt.subplot(111)
for ii in xrange(n_fac):
xx = del_expr_mean_lo[ii]
yy = del_expr_mean_hi[ii]
ax.scatter(xx, yy, lw=0, color=clr[ii], s=10.0)
ax.plot([xx, 1/2*(xx+yy)],
[yy, 1/2*(xx+yy)], color=clr[ii])
min_val, max_val = -3.5, 8
ax.plot([min_val, max_val],
[min_val, max_val], 'k', alpha=0.1)
ax.vlines(0, min_val, max_val, 'r')
ax.hlines(0, min_val, max_val, 'r')
ax.set_xlim([min_val, max_val])
ax.set_ylim([min_val, max_val])
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xlabel('Low Demand: Relative Expression')
ax.set_ylabel('High Demand: Relative Expression')
plt.savefig('{}/{}.Lo_vs_Hi.svg'.format(path_Figures, expr_id))
plt.show()
plt.close()
In [ ]:
%matplotlib inline
coef_ix = np.array(df_cfg['cfg_obs_lut'], dtype=int)
perf_id = 'meanRT'
subgraph_perf = {}
for expr_ii, expr_id in enumerate(df_cfg['cfg_key_label'][()]['Experiment_ID']):
perf_lo = df_perf[expr_id]['lo'][perf_id].mean(axis=-1)
perf_hi = df_perf[expr_id]['hi'][perf_id].mean(axis=-1)
expr_perf = {}
for cnd_ii, cnd_id in enumerate(df_cfg['cfg_key_label'][()]['Condition_ID']):
cnd_perf = {}
for ii, fac_key in enumerate(np.sort(sort_fac_dict.keys())):
fac_ii = sort_fac_dict[fac_key]
sel_fac_coef = df_subgraph['system_subgraph'][fac_ii]['expr_coef'][coef_ix]
pos_expr = sel_fac_coef[:, expr_ii, :, :, :, :][:, cnd_ii, 1, 0, :]
neg_expr = sel_fac_coef[:, expr_ii, :, :, :, :][:, cnd_ii, 1, 1, :]
rel_expr = pos_expr-neg_expr
cnd_perf[fac_key] = stats.spearmanr(rel_expr.mean(axis=-1),
perf_hi-perf_lo)
expr_perf[cnd_id] = cnd_perf
subgraph_perf[expr_id] = expr_perf
In [ ]:
# Subgraph Node Strengths
fac_adj = np.array([df_subgraph['system_subgraph'][sort_fac_dict[fac_key]]['subnet_roi']
for fac_key in np.sort(sort_fac_dict.keys())])
fac_cfg = convert_adj_matr_to_cfg_matr(fac_adj)
subgraph_ns = fac_adj.mean(axis=1)
subgraph_ns_null = []
for pp in xrange(10000):
fac_adj_null = np.array([convert_conn_vec_to_adj_matr(np.random.permutation(cfg))
for cfg in fac_cfg])
subgraph_ns_null.append(fac_adj_null.mean(axis=1))
subgraph_ns_null = np.array(subgraph_ns_null)
# Iterate over Subgraph Performance Predictor
subgraph_partc = {}
for expr_ii, expr_id in enumerate(subgraph_perf.keys()):
expr_partc = {}
for cnd_ii, cnd_id in enumerate(subgraph_perf[expr_id].keys()):
sel_spear = np.array([subgraph_perf[expr_id][cnd_id][fac_key]
for fac_key in np.sort(sort_fac_dict.keys())])
## Performance Plot
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
clr = np.array(['k' for ff in xrange(n_fac)])
for ix, pv in enumerate(sel_spear[:, 1]):
if pv < 0.05:
clr[ix] = 'r'
ax.bar(xrange(n_fac), sel_spear[:, 0], lw=0, color=clr)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticks(np.arange(n_fac)+0.4);
ax.set_xticklabels(np.sort(sort_fac_dict.keys()));
ax.set_xlim([0, n_fac])
ax.set_ylim([-0.55, 0.55])
ax.set_xlabel('Subgraphs')
ax.set_ylabel('rho(Relative Expression, RT Cost)')
plt.savefig('{}/Subgraph_Performance.{}.{}.svg'.format(path_Figures, expr_id, cnd_id))
plt.show()
## Subgraph participation score
expr_partc[cnd_id] = {'real': np.dot(np.arctanh(sel_spear[:, 0]), subgraph_ns),
'null': np.array([np.dot(np.arctanh(sel_spear[:, 0]), s_ns)
for s_ns in subgraph_ns_null])}
subgraph_partc[expr_id] = expr_partc
In [ ]:
from mayavi import mlab
import nibabel as nib
cmap = 'RdBu'
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'][df_to_yeo['sort_laus_to_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))
subgraph_partc_pixmap = {}
for expr_ii, expr_id in enumerate(subgraph_partc.keys()):
xmax = [0.005, 0.01]
for cnd_id in subgraph_partc[expr_id].keys():
real_partc = subgraph_partc[expr_id][cnd_id]['real'].copy()
null_partc = subgraph_partc[expr_id][cnd_id]['null'].copy()
pv = np.mean(np.abs(null_partc) > np.abs(real_partc), axis=0)
is_sig = Echobase.Statistics.FDR.fdr.bhp(pv, alpha=0.05, dependent=True)
real_partc[~is_sig] = 0
# Iterate over hemisphere of the pial surface
for hemi in pial_hemi.keys():
n_vert = len(pial_hemi[hemi]['vert'])
pial_scalars = np.zeros(n_vert)
# Assign subcortical
nonparc_lbl = []
for roi_ii, roi in enumerate(df_to_yeo['df_laus_yeo'][df_to_yeo['sort_laus_to_yeo']]):
lbl_file = label_files[roi_ii]
sys_lbl = roi[3]
if roi[2] != hemi:
continue
if (sys_lbl == 'CRB') or (sys_lbl == 'SUB'):
continue
nonparc_lbl.append([parc_ix for parc_ix in nib.freesurfer.io.read_label(lbl_file)])
nonparc_lbl = np.array(nonparc_lbl)
parc_subcort = np.setdiff1d(np.arange(n_vert), nonparc_lbl)
sys_lbl = df_to_yeo['yeo_lbl'][df_to_yeo['sort_laus_to_yeo']]
pial_scalars[parc_subcort] = np.sum(real_partc[sys_lbl == 'SUB']) / np.sum(real_partc[sys_lbl == 'SUB'] != 0)
# Iterate over brain regions
for roi_ii, roi in enumerate(df_to_yeo['df_laus_yeo'][df_to_yeo['sort_laus_to_yeo']]):
lbl_file = label_files[roi_ii]
sys_lbl = roi[3]
if roi[2] != hemi:
continue
if (sys_lbl == 'CRB') or (sys_lbl == 'SUB'):
continue
# Load the file and add scalar to the vertices
parc_lbl = nib.freesurfer.io.read_label(lbl_file)
pial_scalars[parc_lbl] = real_partc[roi_ii]
# 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=cmap, data_range=[-xmax[expr_ii], xmax[expr_ii]], use_default_range=False)
lut = surf.module_manager.scalar_lut_manager.lut.table.to_array()[::-1, :]
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)
subgraph_partc_pixmap[hemi] = pixmap
fig = plt.figure(figsize=(6,6), dpi=300.0)
ax = fig.add_subplot(2,2,1); ax.imshow(subgraph_partc_pixmap['RH']['Sag_PA']); ax.set_axis_off()
ax = fig.add_subplot(2,2,2); ax.imshow(subgraph_partc_pixmap['LH']['Sag_AP']); ax.set_axis_off()
ax = fig.add_subplot(2,2,3); ax.imshow(subgraph_partc_pixmap['RH']['Sag_AP'][:, ::-1, :]); ax.set_axis_off()
ax = fig.add_subplot(2,2,4); ax.imshow(subgraph_partc_pixmap['LH']['Sag_PA'][:, ::-1, :]); ax.set_axis_off()
fig.savefig('{}/Subgraph_Participation.{}.{}.svg'.format(path_Figures, expr_id, cnd_id))
plt.show()
In [ ]:
plt.figure(figsizedpi=300.0)
ax = plt.subplot(111)
mat = ax.imshow(subgraph_partc_pixmap['RH']['Sag_PA'], cmap='RdBu');
plt.colorbar(mat, ax=ax)
plt.savefig('{}/RdBu_Colorbar.svg'.format(path_Figures))
plt.show()
In [ ]: