Homec4science
Diffusion EEG_dynamic_task_functional_connectivity (master) Browse

EEG_conn/connectivity_group.py
d7c9e3f3df76master

Last Change
Raw File

connectivity_group.py
View Options

import os
import numpy as np
import matplotlib.pyplot as plt
import json
import pandas as pd
import hashlib
import csv
import mne
import scipy.io
from statsmodels.stats.multitest import multipletests
from scipy.stats import mannwhitneyu, shapiro, hmean, pearsonr, spearmanr
from functions import compute_distances, plot_statistical_test_results, \
get_biosemi_to_1020_mapping, compute_connectivity_matrix
param_file = "npm_8-13Hz_0.7s_R.json"
#param_file = "npm_8-13Hz_0.7s_P.json"
#param_file = "npm_14-30Hz_0.7s_R.json"
#param_file = "npm_14-30Hz_0.7s_P.json"
with open(param_file) as json_file:
param_dict = json.load(json_file)
subjects = np.array(sorted(param_dict['subjects'])).astype(str)
n_of_plots_in_line = 3
show_electrodes = 0
do_rename_channels = True
results_base_dir = param_dict['result_dir']
data_base_dir = param_dict['datapath']
reference_raw_file = data_base_dir + "000_P1.vhdr"
conditions = param_dict['conditions']
methods = param_dict['methods']
franges = param_dict['frequency_band']
tmin_tmax = param_dict['tmin_tmax']
csd_params = param_dict['CSD_parameters']
add_label = param_dict['add_label']
min_edge_pvals = []
epochs_tmin = -2
epochs_tmax = 4
reference_raw = mne.io.read_raw_brainvision(reference_raw_file, preload=False)
if do_rename_channels:
names_dict = {}
ch_ind = 0
for prefix in ['A', 'B', 'C', 'D']:
for i in range(0, 32):
name = prefix + str(i + 1)
if name != reference_raw.info['ch_names'][ch_ind]:
names_dict.update({reference_raw.info['ch_names'][ch_ind]: name})
ch_ind += 1
reference_raw.rename_channels(names_dict)
w = csv.writer(open(results_base_dir + "electrodes_dict.csv", "w"))
for key, val in names_dict.items():
w.writerow([key, val])
types_dict = {}
for name in reference_raw.info['ch_names']:
if name[0] != 'E':
types_dict.update({name: 'eeg'})
else:
types_dict.update({name: 'eog'})
reference_raw.set_channel_types(types_dict)
picks = mne.pick_types(reference_raw.info, eeg=True)
n_all_el = picks.shape[0]
my_montage = mne.channels.read_custom_montage("biosemi128.sfp")
reference_raw.set_montage(my_montage)
print("Reference raw loaded")
for condition in conditions:
for f_range in franges:
f_range_low = f_range[0]
f_range_high = f_range[1]
for tmin_tmax_range in tmin_tmax:
if tmin_tmax_range is None:
conn_tmin = -0.2
conn_tmax = 0.5
else:
conn_tmin = tmin_tmax_range[0]
conn_tmax = tmin_tmax_range[1]
for csd_param in csd_params:
if csd_param is None:
do_csd = False
else:
do_csd = True
stiffnes = csd_param[0]
lambd = csd_param[1]
for method in methods:
connectivity_full = np.empty(shape=(n_all_el, n_all_el, len(subjects)))
labels_true = []
subj_pos = 0
n_controls = 0
n_patients = 0
suffix_base = ""
suffix_base += "_" + str(f_range_low) + "-" + str(f_range_high) + "_" + method + "_T_" + \
str(conn_tmin) + "_" + str(conn_tmax) + "s"
if do_csd:
suffix_base += "_csd" + str(stiffnes) + "_" + str(lambd)
else:
suffix_base += "_nocsd"
if condition == "proposer":
suffix_base += "P"
elif condition == "responder":
suffix_base += "R"
for subj in subjects:
if subj[0] == '0':
n_controls += 1
labels_true.append(1)
if condition == "responder":
file_name = subj + "_R1"
elif condition == "proposer":
file_name = subj + "_P1"
else:
raise ValueError("Unrecognised condition" + condition)
elif subj[0] == '1':
n_patients += 1
labels_true.append(2)
if condition == "responder":
file_name = subj + "_R1"
elif condition == "proposer":
file_name = subj + "_P1"
else:
raise ValueError("Unrecognised condition" + condition)
############################# Loading/computing connectivity for each subject #########################
if os.path.exists(results_base_dir + "conn_matrices/" + 'con' + file_name + suffix_base + '.npy'):
print(
"Reading precomputed connectivity matrix:\t" + results_base_dir + "conn_matrices/" + 'con' + file_name + suffix_base + '.npy')
c_matrix = np.load(results_base_dir + "conn_matrices/" + 'con' + file_name + suffix_base + '.npy')
else:
c_matrix = compute_connectivity_matrix(f_range, tmin_tmax_range, csd_param, method,
condition, subj, epochs_tmin, epochs_tmax,
results_base_dir + "conn_matrices/", data_base_dir)
connectivity_full[:, :, subj_pos] = c_matrix
n_channels = c_matrix.shape[0]
subj_pos += 1
labels_true = np.array(labels_true).astype(int)
el_set_index = -1
for selected_channels_param in param_dict['electrode_subset']:
el_set_index += 1
if type(selected_channels_param) is list:
selected_channels = selected_channels_param
elif type(selected_channels_param) is str:
if selected_channels_param == 'All':
selected_channels = ['A1', 'A2', 'A3', 'A4', 'A5', 'A6', 'A7', 'A8', 'A9', 'A10', 'A11',
'A12', 'A13', 'A14', 'A15', 'A16', 'A17', 'A18', 'A19', 'A20', 'A21',
'A22', 'A23', 'A24', 'A25', 'A26', 'A27', 'A28', 'A29', 'A30', 'A31',
'A32', 'B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B9', 'B10',
'B11', 'B12', 'B13', 'B14', 'B15', 'B16', 'B17', 'B18', 'B19', 'B20',
'B21', 'B22', 'B23', 'B24', 'B25', 'B26', 'B27', 'B28', 'B29', 'B30',
'B31', 'B32', 'C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8', 'C9',
'C10', 'C11', 'C12', 'C13', 'C14', 'C15', 'C16', 'C17', 'C18', 'C19',
'C20', 'C21', 'C22', 'C23', 'C24', 'C25', 'C26', 'C27', 'C28', 'C29',
'C30', 'C31', 'C32', 'D1', 'D2', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15', 'D16', 'D17', 'D18',
'D19', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26', 'D27', 'D28',
'D29', 'D30', 'D31', 'D32']
elif selected_channels_param == '32el10-20':
electrode_mapping = get_biosemi_to_1020_mapping(32)
selected_channels = list(electrode_mapping.keys())
elif selected_channels_param == '64el10-20':
electrode_mapping = get_biosemi_to_1020_mapping(64)
selected_channels = list(electrode_mapping.keys())
else:
raise ValueError("Wrong channel selection parameter: should be a list of channel labels or 'All'")
else:
raise ValueError("Wrong channel selection parameter: should be a list of channel labels or 'All'")
if selected_channels_param == '32el10-20' or selected_channels_param == '64el10-20':
el_labels = [electrode_mapping[key] for key in selected_channels]
else:
el_labels = selected_channels
suffix = suffix_base + "_" + str(len(selected_channels)) + "el_"
text_el = "".join(selected_channels)
suffix += hashlib.sha256(text_el.encode()).hexdigest()[:8]
suffix += add_label
electrode_dist_m = compute_distances(selected_channels, ref_data_file=reference_raw_file,
vis=show_electrodes)
n_channels = len(selected_channels)
selected_ch_ind = [mne.pick_channels(reference_raw.info['ch_names'], include=[ch])[0] for ch in selected_channels]
connectivity_subset = connectivity_full[selected_ch_ind, :, :][:, selected_ch_ind, :]
for i in range(connectivity_subset.shape[2]):
connectivity_subset[:, :, i][np.triu_indices_from(connectivity_subset[:, :, i], 1)] = 0
xyz = np.empty((n_channels, 3))
ch_names = np.empty(n_channels)
for i, rawind in enumerate(selected_ch_ind):
xyz[i, :] = reference_raw.info['chs'][rawind]['loc'][0:3]
xyz[:, 1] += -0.017
connectivity_dict = {'conn': connectivity_subset, 'subjects': subjects,
'xyz': xyz, 'el_labels': el_labels,
'design_matrix': np.vstack([labels_true - 1, -1 * (labels_true - 2)])}
scipy.io.savemat(results_base_dir + 'con' + suffix_base + add_label + '.mat', connectivity_dict)
#Non parametric test, gaussianity test
ustats = np.empty_like(connectivity_subset[:, :, 0])
ustats[:] = np.nan
upvalue = np.empty_like(connectivity_subset[:, :, 0])
upvalue[:] = np.nan
isgaussian_stats = np.empty_like(connectivity_subset[:, :, 0])
isgaussian_stats[:] = np.nan
isgaussian_pvalue = np.empty_like(connectivity_subset[:, :, 0])
isgaussian_pvalue[:] = np.nan
flat_mwupvalues = np.empty(shape=(int(n_channels * (n_channels - 1) / 2)))
ind = 0
for i in range(0, n_channels):
for j in range(i + 1, n_channels):
ustats[j, i], upvalue[j, i] = mannwhitneyu(connectivity_subset[j, i, :n_controls], connectivity_subset[j, i, n_controls:])
isgaussian_stats[j, i], isgaussian_pvalue[j, i] = shapiro(connectivity_subset[j, i, :n_controls])
flat_mwupvalues[ind] = upvalue[j, i]
ind += 1
#Multiple compatison correction can only be done on a vector
h_true, corrected_mwupvalues_flat, als, alb = multipletests(flat_mwupvalues, method='fdr_bh')
corrected_mwupvalue = np.empty_like(upvalue)
corrected_mwupvalue[:] = np.nan
ind = 0
for i in range(0, n_channels):
for j in range(i + 1, n_channels):
corrected_mwupvalue[j, i] = corrected_mwupvalues_flat[ind]
ind += 1
df_mwustats = pd.DataFrame(ustats, columns=el_labels, index=el_labels)
df_mwustats.to_csv(results_base_dir + "mwustats" + suffix + '.csv', index=False, header=True, float_format='%.5f')
df_mwupval = pd.DataFrame(upvalue, columns=el_labels, index=el_labels)
df_mwupval.to_csv(results_base_dir + "mwupval" + suffix + '.csv', index=False, header=True, float_format='%.5f')
df_corr_mwupval = pd.DataFrame(corrected_mwupvalue, columns=el_labels, index=el_labels)
df_corr_mwupval.to_csv(results_base_dir + "mwupval_corrected_" + suffix + '.csv', index=False,
header=True, float_format='%.5f')
min_edge_pvals.append([len(selected_channels), suffix[-8:], method, f_range[0], f_range[1],
conn_tmin, conn_tmax, stiffnes, lambd, -1 * np.log10(np.nanmin(corrected_mwupvalue))])
mean_pval = np.zeros(upvalue.shape[0])
symm_upvalue = np.copy(upvalue)
symm_upvalue = np.nan_to_num(symm_upvalue)
symm_upvalue += symm_upvalue.T
for i in range(corrected_mwupvalue.shape[0]):
if i == 0:
mean_pval[i] = hmean(symm_upvalue[i, i + 1:])
elif 0 < i < symm_upvalue.shape[0] - 1:
mean_pval[i] = hmean(np.hstack([symm_upvalue[i, 0:i], symm_upvalue[i, i + 1:]]))
else:
mean_pval[i] = hmean(symm_upvalue[i, 0:i])
df_mean_mwupval = pd.DataFrame(mean_pval, index=selected_channels)
df_mean_mwupval.to_csv(results_base_dir + "mean_u_pvalues" + suffix + '.csv',
index=selected_channels, header=["mean_pval"], float_format='%.5f')
plot_statistical_test_results(isgaussian_stats, isgaussian_pvalue, labels=el_labels,
writepath=results_base_dir + "gaussianity_test" + suffix + ".png",
title='Gaussianity test, controls, ' + method)
plot_statistical_test_results(ustats, upvalue, labels=el_labels,
writepath=results_base_dir + "mannwhitney" + suffix + ".png",
title='Mann-Whitney rank test, ' + method, left_title="Mann-Whitney U statistic")
plot_statistical_test_results(ustats, corrected_mwupvalue, labels=el_labels,
writepath=results_base_dir + "mannwhitney_corr" + suffix + ".png",
title='Mann-Whitney rank test, corrected, ' + method, left_title="Mann-Whitney U statistic")
plot_statistical_test_results(ustats, -1*np.log10(corrected_mwupvalue), labels=el_labels,
writepath=results_base_dir + "mannwhitney_corr_log" + suffix + ".png",
title='-log Mann-Whitney rank test, corrected, ' + method,
left_title="Mann-Whitney U statistic")
results_dir = results_base_dir + suffix[1:] + "/"
if not os.path.exists(results_dir):
os.mkdir(results_dir)
print("Results saved at: " + results_dir)
el_dist_flat = electrode_dist_m[np.tril_indices(connectivity_subset.shape[0], -1)]
############## plotting distance matrix ############################
fig = plt.figure()
im = plt.imshow(electrode_dist_m, cmap='hot')
step = max(1, int(np.log2(len(selected_channels)+1) / 2))
plt.xticks(np.arange(len(selected_channels), step=step), selected_channels[::step], rotation=270)
plt.yticks(np.arange(len(selected_channels), step=step), selected_channels[::step])
plt.colorbar(im, fraction=0.046, pad=0.04)
plt.savefig(results_dir + "distance_matrix" + suffix + ".png", dpi=600)
plt.close(fig)
result_flat = np.zeros(shape=(len(subjects), len(np.tril_indices(connectivity_subset.shape[0], -1)[0])))
for i in range(0, len(subjects)):
result_flat[i, :] = connectivity_subset[:, :, i][
np.tril_indices(connectivity_subset.shape[0], -1)]
for i in range(connectivity_subset.shape[2]):
connectivity_subset[:, :, i][np.triu_indices_from(connectivity_subset[:, :, i], 1)] = np.nan
for p_number in range(0, int(max(n_patients, n_controls) / n_of_plots_in_line)):
fig = plt.figure('Connectivity ' + str(p_number) + ', ' + method, figsize=[12, 6])
for i in range(0, n_of_plots_in_line):
ax = plt.subplot(2, n_of_plots_in_line, i + 1)
if p_number * n_of_plots_in_line + i < n_controls:
plt.title(subjects[p_number * n_of_plots_in_line + i])
im = plt.imshow(connectivity_subset[:, :, p_number * n_of_plots_in_line + i], cmap='twilight_shifted')
plt.xticks(np.arange(len(el_labels), step=step), el_labels[::step], rotation=270)
plt.yticks(np.arange(len(el_labels), step=step), el_labels[::step])
for i in range(0, n_of_plots_in_line):
ax = plt.subplot(2, n_of_plots_in_line, i + 1 + n_of_plots_in_line)
if p_number * n_of_plots_in_line + i < n_patients:
plt.title(subjects[n_controls + p_number * n_of_plots_in_line + i])
im = plt.imshow(connectivity_subset[:, :, n_controls + p_number * n_of_plots_in_line + i], cmap='twilight_shifted')
plt.xticks(np.arange(len(el_labels), step=step), el_labels[::step], rotation=270)
plt.yticks(np.arange(len(el_labels), step=step), el_labels[::step])
plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)
cax = plt.axes([0.85, 0.12, 0.03, 0.75])
cbar = plt.colorbar(cax=cax)
cbar.ax.set_ylabel('Connectivity strength')
plt.subplots_adjust(top=0.935, bottom=0.075, left=0.04, right=0.835, hspace=0.3, wspace=0.254)
plt.savefig(results_dir + "connectiviy" + str(p_number) + suffix + ".png", dpi=600)
plt.close(fig)
fig = plt.figure(
'Connectivity (y-axis) vs distance (x-axis) correlation ' + str(p_number) + ', ' + method, figsize=[12, 6])
for i in range(0, n_of_plots_in_line):
ax = plt.subplot(2, n_of_plots_in_line, i + 1)
if p_number * n_of_plots_in_line + i < n_controls:
corr, pcorr = spearmanr(100*el_dist_flat, result_flat[p_number * n_of_plots_in_line + i, :])
plt.title(subjects[p_number * n_of_plots_in_line + i] + "\ncorr=" +
format(corr, '.2f') + " p=" + format(pcorr, '.2f'))
plt.scatter(100*el_dist_flat, result_flat[p_number * n_of_plots_in_line + i, :])
#plt.xlabel("Distance (cm)")
plt.ylabel("ImCoh")
axes = plt.gca()
for i in range(0, n_of_plots_in_line):
ax = plt.subplot(2, n_of_plots_in_line, i + 1 + n_of_plots_in_line)
if p_number * n_of_plots_in_line + i < n_patients:
corr, pcorr = pearsonr(100 * el_dist_flat,
result_flat[p_number * n_of_plots_in_line + i, :])
plt.title(subjects[n_controls+p_number * n_of_plots_in_line + i] +
"\ncorr=" + format(corr, '.2f') + ", p=" + format(pcorr, '.2f'))
plt.scatter(100*el_dist_flat, result_flat[n_controls + p_number * n_of_plots_in_line + i, :])
plt.xlabel("Distance (cm)")
plt.ylabel("ImCoh")
axes = plt.gca()
plt.tight_layout(pad=1)
plt.savefig(results_dir + "conn_vs_dist" + str(p_number) + suffix + ".png", dpi=600)
plt.close(fig)
for p_number in range(0, int(max(n_patients, n_controls) / n_of_plots_in_line)):
fig = plt.figure('Weight distribution, ' + str(p_number) + ', ' + method, figsize=[12, 6])
for i in range(0, n_of_plots_in_line):
ax = plt.subplot(2, n_of_plots_in_line, i + 1)
if p_number * n_of_plots_in_line + i < n_controls:
plt.title(subjects[p_number * n_of_plots_in_line + i])
im = plt.hist(connectivity_subset[:, :, p_number * n_of_plots_in_line + i][np.where(
connectivity_subset[:, :, p_number * n_of_plots_in_line + i] != 0)].flatten())
axes = plt.gca()
axes.set_ylim([0, 100])
axes.set_xlim([np.nanmin(connectivity_subset), np.nanmax(connectivity_subset)])
for i in range(0, n_of_plots_in_line):
ax = plt.subplot(2, n_of_plots_in_line, i + 1 + n_of_plots_in_line)
if p_number * n_of_plots_in_line + i < n_patients:
plt.title(subjects[n_controls + p_number * n_of_plots_in_line + i])
im = plt.hist(connectivity_subset[:, :, n_controls + p_number * n_of_plots_in_line + i][np.where(
connectivity_subset[:, :, n_controls + p_number * n_of_plots_in_line + i] != 0)].flatten())
axes = plt.gca()
axes.set_ylim([0, 100])
axes.set_xlim([np.nanmin(connectivity_subset), np.nanmax(connectivity_subset)])
plt.tight_layout(pad=1)
plt.savefig(results_dir + "weight_distr" + str(p_number) + suffix + ".png", dpi=600)
plt.close(fig)
min_edge_pvals_df = pd.DataFrame(min_edge_pvals, columns=['n-el', 'el-hash', 'method', 'freq-min', 'freq-max',
"tmin", "tmax", "csd-st", "csd-lam", '-log10(min-pval)'])
np.savetxt(results_base_dir + param_file[:-5] + 'min_edge_pvals.csv', min_edge_pvals_df.values,
fmt=['%d', '%s', '%s', '%d', '%d', '%.2f', '%.2f', '%.2f', '%.2e', '%.2f'], delimiter=', ',
header=''+','.join(min_edge_pvals_df.columns), comments='')
c4science · Help