In [ ]:
try:
%load_ext autoreload
%autoreload 2
%reset
except:
print 'NOT IPYTHON'
from __future__ import division
import os
import sys
import glob
import numpy as np
import pandas as pd
import seaborn as sns
import scipy.stats as stats
import statsmodels.api as sm
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_ExpData = path_PeriphData + '/e01-FuncNetw'
path_AtlasData = path_Remotes + '/CORE.MRI_Atlases'
path_Figures = './e01-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 BOLD
df_navon = io.loadmat('{}/NavonBlockedSeriesScale125.mat'.format(path_CoreData), struct_as_record=False)
df_stroop = io.loadmat('{}/StroopBlockedSeriesScale125.mat'.format(path_CoreData), struct_as_record=False)
n_subj = 28
n_fix_block = 12 # Disregard the final fixation block
n_tsk_block = 6
n_roi = 262
bad_roi = [242]
n_good_roi = n_roi-len(bad_roi)
# 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]}
# 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)
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, :]}
}
}
In [ ]:
def comp_fconn(bold, alpha=0.05, dependent=False):
n_roi, n_tr = bold.shape
adj = np.arctanh(np.corrcoef(bold))
cfg_vec = convert_adj_matr_to_cfg_matr(adj.reshape(-1, n_roi, n_roi))[0, :]
# Separate edges based on sign
cfg_vec_pos = cfg_vec.copy()
cfg_vec_pos[cfg_vec_pos < 0] = 0
cfg_vec_neg = -1*cfg_vec.copy()
cfg_vec_neg[cfg_vec_neg < 0] = 0
adj_pos = convert_conn_vec_to_adj_matr(cfg_vec_pos)
adj_neg = convert_conn_vec_to_adj_matr(cfg_vec_neg)
return adj_pos, adj_neg
In [ ]:
for subj_id in xrange(n_subj):
proc_item = '{}/Subject_{}.Navon'.format(path_ExpData, subj_id)
print(proc_item)
adj_dict = {'lo': {'fix': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))},
'task': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))}
},
'hi': {'fix': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))},
'task': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))},
}}
# Process Fixation Blocks
cnt = 0
for fix_block in xrange(n_fix_block):
data = np.array(df_navon['data'][subj_id][fix_block].NFix, dtype='f').T
data = data[np.setdiff1d(np.arange(n_roi), bad_roi), :]
if (fix_block % 2) == 0:
adj_dict['lo']['fix']['pos'][cnt, :, :], adj_dict['lo']['fix']['neg'][cnt, :, :] = comp_fconn(data)
if (fix_block % 2) == 1:
adj_dict['hi']['fix']['pos'][cnt, :, :], adj_dict['hi']['fix']['neg'][cnt, :, :] = comp_fconn(data)
cnt += 1
# Process Task Blocks
cnt = 0
for tsk_block in xrange(n_tsk_block):
# Low demand
data = np.array(df_navon['data'][subj_id][tsk_block].NS, dtype='f').T
data = data[np.setdiff1d(np.arange(n_roi), bad_roi), :]
adj_dict['lo']['task']['pos'][cnt, :, :], adj_dict['lo']['task']['neg'][cnt, :, :] = comp_fconn(data)
# High demand
data = np.array(df_navon['data'][subj_id][tsk_block].S, dtype='f').T
data = data[np.setdiff1d(np.arange(n_roi), bad_roi), :]
adj_dict['hi']['task']['pos'][cnt, :, :], adj_dict['hi']['task']['neg'][cnt, :, :] = comp_fconn(data)
cnt += 1
np.savez(proc_item, adj_dict=adj_dict)
In [ ]:
for subj_id in xrange(n_subj):
proc_item = '{}/Subject_{}.Stroop'.format(path_ExpData, subj_id)
print(proc_item)
adj_dict = {'lo': {'fix': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))},
'task': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))}
},
'hi': {'fix': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))},
'task': {'pos': np.zeros((n_tsk_block, n_good_roi, n_good_roi)),
'neg': np.zeros((n_tsk_block, n_good_roi, n_good_roi))},
}}
# Process Fixation Blocks
cnt = 0
for fix_block in xrange(n_fix_block):
data = np.array(df_stroop['data'][subj_id][fix_block].SFix, dtype='f').T
data = data[np.setdiff1d(np.arange(n_roi), bad_roi), :]
if (fix_block % 2) == 0:
adj_dict['lo']['fix']['pos'][cnt, :, :], adj_dict['lo']['fix']['neg'][cnt, :, :] = comp_fconn(data)
if (fix_block % 2) == 1:
adj_dict['hi']['fix']['pos'][cnt, :, :], adj_dict['hi']['fix']['neg'][cnt, :, :] = comp_fconn(data)
cnt += 1
# Process Task Blocks
cnt = 0
for tsk_block in xrange(n_tsk_block):
# Low demand
data = np.array(df_stroop['data'][subj_id][tsk_block].IE, dtype='f').T
data = data[np.setdiff1d(np.arange(n_roi), bad_roi), :]
adj_dict['lo']['task']['pos'][cnt, :, :], adj_dict['lo']['task']['neg'][cnt, :, :] = comp_fconn(data)
# High demand
data = np.array(df_stroop['data'][subj_id][tsk_block].E, dtype='f').T
data = data[np.setdiff1d(np.arange(n_roi), bad_roi), :]
adj_dict['hi']['task']['pos'][cnt, :, :], adj_dict['hi']['task']['neg'][cnt, :, :] = comp_fconn(data)
cnt += 1
np.savez(proc_item, adj_dict=adj_dict)
In [ ]:
expr_dict = {}
for expr_id in ['Stroop', 'Navon']:
df_list = glob.glob('{}/Subject_*.{}.npz'.format(path_ExpData, expr_id))
for df_subj in df_list:
subj_id = int(df_subj.split('/')[-1].split('.')[0].split('_')[1])
if subj_id not in expr_dict.keys():
expr_dict[subj_id] = {}
expr_dict[subj_id][expr_id] = df_subj
In [ ]:
# Generate a dictionary of all key names
cfg_key_names = ['Subject_ID', 'Experiment_ID', 'Condition_ID', 'Task_ID', 'CorSign_ID', 'Block_ID']
cfg_key_label = {'Subject_ID': np.arange(n_subj),
'Experiment_ID': ['Stroop', 'Navon'],
'Condition_ID': ['lo', 'hi'],
'Task_ID': ['fix', 'task'],
'CorSign_ID': ['pos', 'neg'],
'Block_ID': np.arange(n_tsk_block)}
cfg_obs_lut = np.zeros((len(cfg_key_label[cfg_key_names[0]]),
len(cfg_key_label[cfg_key_names[1]]),
len(cfg_key_label[cfg_key_names[2]]),
len(cfg_key_label[cfg_key_names[3]]),
len(cfg_key_label[cfg_key_names[4]]),
len(cfg_key_label[cfg_key_names[5]])))
# Iterate over all cfg key labels and generate a LUT matrix and a config matrix
key_cnt = 0
cfg_matr = []
for key_0_ii, key_0_id in enumerate(cfg_key_label[cfg_key_names[0]]):
for key_1_ii, key_1_id in enumerate(cfg_key_label[cfg_key_names[1]]):
adj_dict = np.load(expr_dict[key_0_id][key_1_id])['adj_dict'][()]
for key_2_ii, key_2_id in enumerate(cfg_key_label[cfg_key_names[2]]):
for key_3_ii, key_3_id in enumerate(cfg_key_label[cfg_key_names[3]]):
for key_4_ii, key_4_id in enumerate(cfg_key_label[cfg_key_names[4]]):
for key_5_ii, cfg_vec in enumerate(convert_adj_matr_to_cfg_matr(adj_dict[key_2_id][key_3_id][key_4_id])):
cfg_obs_lut[key_0_ii, key_1_ii, key_2_ii,
key_3_ii, key_4_ii, key_5_ii] = key_cnt
cfg_matr.append(cfg_vec)
key_cnt += 1
cfg_matr = np.array(cfg_matr)
cfg_matr_orig = cfg_matr.copy()
# Normalize sum of edge weights to 1
cfg_L1 = np.linalg.norm(cfg_matr, axis=1, ord=1)
cfg_L1[cfg_L1 == 0] = 1.0
cfg_matr = (cfg_matr.T / cfg_L1).T
# Rescale edge weight to unit L2-Norm
cfg_L2 = np.zeros_like(cfg_matr)
for subj_ii in xrange(len(cfg_key_label['Subject_ID'])):
grp_ix = np.array(cfg_obs_lut[subj_ii, :, :, :, :, :].reshape(-1), dtype=int)
cfg_L2[grp_ix, :] = np.linalg.norm(cfg_matr[grp_ix, :], axis=0, ord=2)
cfg_L2[cfg_L2 == 0] = 1.0
cfg_matr = cfg_matr / cfg_L2
np.savez('{}/Population.Configuration_Matrix.npz'.format(path_ExpData),
cfg_matr_orig=cfg_matr_orig,
cfg_matr=cfg_matr,
cfg_L2=cfg_L2,
cfg_obs_lut=cfg_obs_lut,
cfg_key_label=cfg_key_label,
cfg_key_names=cfg_key_names)
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr_orig']
n_grp = len(df['cfg_key_label'][()]['Subject_ID'])
grp_edge_wt = []
for grp_ii in xrange(n_grp):
grp_ix = np.array(cfg_obs_lut[grp_ii, :, :, :, :, :].reshape(-1), dtype=int)
grp_edge_wt.append(np.mean(cfg_matr[grp_ix, :], axis=1))
grp_edge_wt = np.array(grp_edge_wt)
mean_grp_edge_wt = np.mean(grp_edge_wt, axis=1)
grp_ord_ix = np.argsort(mean_grp_edge_wt)[::-1]
### Plot Subject Distribution
print(stats.f_oneway(*(grp_edge_wt)))
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
bp = ax.boxplot(grp_edge_wt[grp_ord_ix, :].T, sym='', patch_artist=True)
Echobase.Plotting.fig_format.set_box_color(bp, [0.0, 0.0, 0.0], [[0.2, 0.2, 0.2] for iii in xrange(n_grp)])
ax.set_ylim(ymin=0)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticklabels([])
ax.set_xlabel('Subjects')
ax.set_ylabel('Weighted Edge Density')
plt.savefig('{}/Wgt_Edge_Density.Subjects.svg'.format(path_Figures))
plt.show()
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr_orig']
n_grp = len(df['cfg_key_label'][()]['CorSign_ID'])
n_subj = len(df['cfg_key_label'][()]['Subject_ID'])
grp_edge_wt = np.zeros((n_grp, n_subj))
for grp_ii in xrange(n_grp):
for subj_ii in xrange(n_subj):
grp_ix = np.array(cfg_obs_lut[subj_ii, :, :, :, :, :][:, :, :, grp_ii, :].reshape(-1), dtype=int)
grp_edge_wt[grp_ii, subj_ii] = np.mean(np.mean(cfg_matr[grp_ix, :], axis=1))
print(stats.ttest_rel(*(grp_edge_wt)))
""
### Plot
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
bp = ax.boxplot(grp_edge_wt.T, patch_artist=True)
Echobase.Plotting.fig_format.set_box_color(bp, [0.0, 0.0, 0.0], [[0.2, 0.2, 0.2] for iii in xrange(n_grp)])
ax.set_ylim(ymin=0)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticklabels(df['cfg_key_label'][()]['CorSign_ID'])
ax.set_xlabel('')
ax.set_ylabel('Weighted Edge Density')
plt.savefig('{}/Wgt_Edge_Density.CorSign.svg'.format(path_Figures))
plt.show()
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr_orig']
n_grp = len(df['cfg_key_label'][()]['Task_ID'])
n_subj = len(df['cfg_key_label'][()]['Subject_ID'])
grp_edge_wt = np.zeros((n_grp, n_subj))
for grp_ii in xrange(n_grp):
for subj_ii in xrange(n_subj):
grp_ix = np.array(cfg_obs_lut[subj_ii, :, :, :, :, :][:, :, grp_ii, :, :].reshape(-1), dtype=int)
grp_edge_wt[grp_ii, subj_ii] = np.mean(np.mean(cfg_matr[grp_ix, :], axis=1))
print(stats.ttest_rel(*(grp_edge_wt)))
### Plot
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
bp = ax.boxplot(grp_edge_wt.T, patch_artist=True)
Echobase.Plotting.fig_format.set_box_color(bp, [0.0, 0.0, 0.0], [[0.2, 0.2, 0.2] for iii in xrange(n_grp)])
ax.set_ylim(ymin=0)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticklabels(df['cfg_key_label'][()]['Task_ID'])
ax.set_xlabel('')
ax.set_ylabel('Weighted Edge Density')
plt.savefig('{}/Wgt_Edge_Density.Task.svg'.format(path_Figures))
plt.show()
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr_orig']
n_grp = len(df['cfg_key_label'][()]['Condition_ID'])
n_subj = len(df['cfg_key_label'][()]['Subject_ID'])
grp_edge_wt = np.zeros((n_grp, n_subj))
for grp_ii in xrange(n_grp):
for subj_ii in xrange(n_subj):
grp_ix = np.array(cfg_obs_lut[subj_ii, :, :, :, :, :][:, grp_ii, :, :, :].reshape(-1), dtype=int)
grp_edge_wt[grp_ii, subj_ii] = np.mean(np.mean(cfg_matr[grp_ix, :], axis=1))
print(stats.ttest_rel(*(grp_edge_wt)))
### Plot
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
bp = ax.boxplot(grp_edge_wt.T, patch_artist=True)
Echobase.Plotting.fig_format.set_box_color(bp, [0.0, 0.0, 0.0], [[0.2, 0.2, 0.2] for iii in xrange(n_grp)])
ax.set_ylim(ymin=0)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticklabels(df['cfg_key_label'][()]['Condition_ID'])
ax.set_xlabel('')
ax.set_ylabel('Weighted Edge Density')
plt.savefig('{}/Wgt_Edge_Density.Condition.svg'.format(path_Figures))
plt.show()
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr_orig']
n_grp = len(df['cfg_key_label'][()]['Experiment_ID'])
n_subj = len(df['cfg_key_label'][()]['Subject_ID'])
grp_edge_wt = np.zeros((n_grp, n_subj))
for grp_ii in xrange(n_grp):
for subj_ii in xrange(n_subj):
grp_ix = np.array(cfg_obs_lut[subj_ii, :, :, :, :, :][grp_ii, :, :, :, :].reshape(-1), dtype=int)
grp_edge_wt[grp_ii, subj_ii] = np.mean(np.mean(cfg_matr[grp_ix, :], axis=1))
print(stats.ttest_rel(*(grp_edge_wt)))
### Plot
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
bp = ax.boxplot(grp_edge_wt.T, patch_artist=True)
Echobase.Plotting.fig_format.set_box_color(bp, [0.0, 0.0, 0.0], [[0.2, 0.2, 0.2] for iii in xrange(n_grp)])
ax.set_ylim(ymin=0)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticklabels(df['cfg_key_label'][()]['Experiment_ID'])
ax.set_xlabel('')
ax.set_ylabel('Weighted Edge Density')
plt.savefig('{}/Wgt_Edge_Density.Expriment.svg'.format(path_Figures))
plt.show()
In [ ]:
perf_stroop_hi = df_perf['Stroop']['hi']['meanRT'].mean(axis=1)
perf_stroop_lo = df_perf['Stroop']['lo']['meanRT'].mean(axis=1)
perf_stroop_cost = perf_stroop_hi-perf_stroop_lo
perf_navon_hi = df_perf['Navon']['hi']['meanRT'].mean(axis=1)
perf_navon_lo = df_perf['Navon']['lo']['meanRT'].mean(axis=1)
perf_navon_cost = perf_navon_hi-perf_navon_lo
print(stats.ttest_rel(perf_stroop_cost, perf_navon_cost))
### Plot
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
bp = ax.boxplot([perf_stroop_cost, perf_navon_cost], patch_artist=True)
Echobase.Plotting.fig_format.set_box_color(bp, [0.0, 0.0, 0.0], [[0.2, 0.2, 0.2] for iii in xrange(2)])
#ax.set_ylim(ymin=0)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
ax.set_xticklabels(['Stroop', 'Navon'])
ax.set_xlabel('')
ax.set_ylabel('Reaction Time Cost (Hi-Lo)')
plt.savefig('{}/RT_Cost.Expriment.svg'.format(path_Figures))
plt.show()
In [ ]:
import nibabel as nib
df_yeo_atlas = nib.load('{}/Yeo_JNeurophysiol11_MNI152/Yeo2011_7Networks_MNI152_FreeSurferConformed1mm_LiberalMask.nii.gz'.format(path_AtlasData))
yeo_matr = df_yeo_atlas.get_data()[..., 0]
yeo_roi = np.unique(yeo_matr)[1:]
yeo_names = ['VIS', 'SMN', 'DAN', 'VAN', 'LIM', 'FPN', 'DMN']
yeo_xyz = {}
M = df_yeo_atlas.affine[:3, :3]
abc = df_yeo_atlas.affine[:3, 3]
for yeo_id in yeo_roi:
yeo_ijk = np.array(np.nonzero(yeo_matr == yeo_id)).T
yeo_xyz[yeo_id] = M.dot(yeo_ijk.T).T + abc.T
df_laus_atlas = nib.load('{}/Lausanne/ROIv_scale125_dilated.nii.gz'.format(path_AtlasData))
laus_matr = df_laus_atlas.get_data()
laus_roi = np.unique(laus_matr)[1:]
laus_xyz = {}
M = df_laus_atlas.affine[:3, :3]
abc = df_laus_atlas.affine[:3, 3]
for laus_id in laus_roi:
laus_ijk = np.array(np.nonzero(laus_matr == laus_id)).T
laus_xyz[laus_id] = M.dot(laus_ijk.T).T + abc.T
laus_yeo_assign = []
for laus_id in laus_roi:
dists = []
for yeo_id in yeo_roi:
dists.append(np.min(np.sum((yeo_xyz[yeo_id] - laus_xyz[laus_id].mean(axis=0))**2, axis=1)))
laus_yeo_assign.append(yeo_names[np.argmin(dists)])
laus_yeo_assign = np.array(laus_yeo_assign)
pd.DataFrame(laus_yeo_assign).to_csv('{}/Lausanne/ROIv_scale125_dilated.Yeo2011_7Networks_MNI152.csv'.format(path_AtlasData))
# Manually replaced subcortical and cerebellar structures as SUB and CBR, respectively.
In [ ]:
# Read in Yeo Atlas
df_laus_yeo = pd.read_csv('{}/LausanneScale125.csv'.format(path_CoreData))
df_laus_yeo = df_laus_yeo[df_laus_yeo.Label_ID != bad_roi[0]+1]
system_lbl = np.array(df_laus_yeo['Yeo2011_7Networks'].as_matrix())
system_name = np.unique(df_laus_yeo['Yeo2011_7Networks'])
n_system = len(system_name)
n_roi = len(system_lbl)
triu_ix, triu_iy = np.triu_indices(n_roi, k=1)
sys_triu_ix, sys_triu_iy = np.triu_indices(n_system, k=0)
# Reorder System Labels and Count ROIs per System
system_srt_ix = np.argsort(system_lbl)
system_cnt = np.array([len(np.flatnonzero(system_lbl == sys_name))
for sys_name in system_name])
system_demarc = np.concatenate(([0], np.cumsum(system_cnt)))
np.savez('{}/Lausanne125_to_Yeo.npz'.format(path_ExpData),
df_laus_yeo=df_laus_yeo,
yeo_lbl=system_lbl,
yeo_name=system_name,
sort_laus_to_yeo=system_srt_ix,
yeo_adj_demarc=system_demarc,
laus_triu=np.triu_indices(n_roi, k=1),
yeo_triu=np.triu_indices(n_system, k=0))
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr']
df_to_yeo = np.load('{}/Lausanne125_to_Yeo.npz'.format(path_ExpData))
n_laus = len(df_to_yeo['yeo_lbl'])
plt.figure(figsize=(5,5));
cnt = 0
for expr_ii, expr_id in enumerate(df['cfg_key_label'][()]['Experiment_ID']):
for sgn_ii, sgn_id in enumerate(df['cfg_key_label'][()]['CorSign_ID']):
grp_ix = np.array(cfg_obs_lut[:, expr_ii, :, :, :, :][:, :, :, sgn_ii, :].reshape(-1), dtype=int)
sel_cfg_matr = cfg_matr[grp_ix, :].mean(axis=0)
adj = convert_conn_vec_to_adj_matr(sel_cfg_matr)
adj_yeo = adj[df_to_yeo['sort_laus_to_yeo'], :][:, df_to_yeo['sort_laus_to_yeo']]
# Plot
ax = plt.subplot(2, 2, cnt+1)
mat = ax.matshow(adj_yeo,
cmap='magma', vmin=0.025)
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.5)
ax.hlines(xx, 0, n_laus, color='w', lw=0.5)
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=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)
ax.set_title('{}-{}'.format(expr_id, sgn_id), fontsize=5.0)
cnt += 1
plt.show()
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_matr = df['cfg_matr']
# Compute Brain System Adjacency Matrices
sys_adj_matr = np.zeros((cfg_matr.shape[0], n_system, n_system))
for sys_ii, (sys_ix, sys_iy) in enumerate(zip(sys_triu_ix, sys_triu_iy)):
sys1 = system_name[sys_ix]
sys2 = system_name[sys_iy]
sys1_ix = np.flatnonzero(system_lbl[triu_ix] == sys1)
sys2_iy = np.flatnonzero(system_lbl[triu_iy] == sys2)
inter_sys_ii = np.intersect1d(sys1_ix, sys2_iy)
if len(inter_sys_ii) == 0:
sys1_ix = np.flatnonzero(system_lbl[triu_ix] == sys2)
sys2_iy = np.flatnonzero(system_lbl[triu_iy] == sys1)
inter_sys_ii = np.intersect1d(sys1_ix, sys2_iy)
mean_conn_sys1_sys2 = np.mean(cfg_matr[:, inter_sys_ii], axis=1)
sys_adj_matr[:, sys_ix, sys_iy] = mean_conn_sys1_sys2
sys_adj_matr[:, sys_iy, sys_ix] = mean_conn_sys1_sys2
np.savez('{}/Full_Adj.Yeo2011_7Networks.npz'.format(path_ExpData),
sys_adj_matr=sys_adj_matr,
cfg_obs_lut=df['cfg_obs_lut'],
cfg_key_label=df['cfg_key_label'],
cfg_key_names=df['cfg_key_names'])
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr']
df_to_yeo = np.load('{}/Lausanne125_to_Yeo.npz'.format(path_ExpData))
n_laus = len(df_to_yeo['yeo_lbl'])
coef_ix = np.array(cfg_obs_lut, dtype=int)
cfg_matr_reshape = cfg_matr[coef_ix, :]
for sgn_ii, sgn_id in enumerate(df['cfg_key_label'][()]['CorSign_ID']):
sel_cfg_matr = (cfg_matr_reshape[:, :, :, 1, sgn_ii, :, :]).mean(axis=-2).mean(axis=-2)
sel_cfg_matr_tv = np.nan*np.zeros(cfg_matr.shape[1])
sel_cfg_matr_pv = np.nan*np.zeros(cfg_matr.shape[1])
for cc in xrange(cfg_matr.shape[1]):
tv, pv = stats.ttest_rel(*sel_cfg_matr[:, :, cc].T)
mean_stroop = np.mean(sel_cfg_matr[:, :, cc], axis=0)[0]
mean_navon = np.mean(sel_cfg_matr[:, :, cc], axis=0)[1]
dv = (mean_stroop - mean_navon) / np.std(sel_cfg_matr[:, :, cc].reshape(-1))
sel_cfg_matr_tv[cc] = dv
sel_cfg_matr_pv[cc] = pv
sig_pv = Echobase.Statistics.FDR.fdr.bhp(sel_cfg_matr_pv, alpha=0.05, dependent=True)
sel_cfg_matr_tv[sig_pv == False] = 0.0
adj = convert_conn_vec_to_adj_matr(sel_cfg_matr_tv)
adj_yeo = adj[df_to_yeo['sort_laus_to_yeo'], :][:, df_to_yeo['sort_laus_to_yeo']]
adj_yeo[np.diag_indices_from(adj_yeo)] = np.nan
# Plot
plt.figure(figsize=(3,3), dpi=300.0)
ax = plt.subplot(111)
mat = ax.matshow(adj_yeo,
cmap='PuOr', vmin=-1.0, vmax=1.0)
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='k', lw=0.5)
ax.hlines(xx, 0, n_laus, color='k', lw=0.5)
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=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)
plt.savefig('{}/Contrast.Expr.{}.svg'.format(path_Figures, sgn_id))
plt.show()
In [ ]:
df = np.load('{}/Population.Configuration_Matrix.npz'.format(path_ExpData))
cfg_obs_lut = df['cfg_obs_lut']
cfg_matr = df['cfg_matr']
df_to_yeo = np.load('{}/Lausanne125_to_Yeo.npz'.format(path_ExpData))
n_laus = len(df_to_yeo['yeo_lbl'])
for expr_ii, expr_id in enumerate(df['cfg_key_label'][()]['Experiment_ID']):
for sgn_ii, sgn_id in enumerate(df['cfg_key_label'][()]['CorSign_ID']):
coef_ix = np.array(cfg_obs_lut, dtype=int)
cfg_matr_reshape = cfg_matr[coef_ix, :]
sel_cfg_matr = cfg_matr_reshape[:, expr_ii, :, 1, sgn_ii, :, :].mean(axis=-2)
sel_cfg_matr_tv = np.nan*np.zeros(cfg_matr.shape[1])
sel_cfg_matr_pv = np.nan*np.zeros(cfg_matr.shape[1])
for cc in xrange(cfg_matr.shape[1]):
tv, pv = stats.ttest_rel(*sel_cfg_matr[:, :, cc].T)
mean_lo = np.mean(sel_cfg_matr[:, :, cc], axis=0)[0]
mean_hi = np.mean(sel_cfg_matr[:, :, cc], axis=0)[1]
dv = (mean_hi - mean_lo) / np.std(sel_cfg_matr[:, :, cc].reshape(-1))
sel_cfg_matr_tv[cc] = dv
sel_cfg_matr_pv[cc] = pv
sig_pv = Echobase.Statistics.FDR.fdr.bhp(sel_cfg_matr_pv, alpha=0.05, dependent=True)
sel_cfg_matr_tv[sig_pv == False] = np.nan
adj = convert_conn_vec_to_adj_matr(sel_cfg_matr_tv)
adj_yeo = adj[df_to_yeo['sort_laus_to_yeo'], :][:, df_to_yeo['sort_laus_to_yeo']]
adj_yeo[np.diag_indices_from(adj_yeo)] = np.nan
# Plot
plt.figure(figsize=(3,3), dpi=300)
ax = plt.subplot(111)
mat = ax.matshow(adj_yeo,
cmap='coolwarm', vmin=-0.5, vmax=0.5)
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='k', lw=0.5)
ax.hlines(xx, 0, n_laus, color='k', lw=0.5)
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=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)
ax.set_title('{}-{}'.format(expr_id, sgn_id), fontsize=5.0)
plt.savefig('{}/Contrast.{}.{}.Hi_Lo.svg'.format(path_Figures, expr_id, sgn_id))
plt.show()
In [ ]: