123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810 |
- import os
- import pandas as pd
- import numpy as np
- import matplotlib.pylab as plt
- import seaborn as sns
- import mr_utilities
- from scipy import stats
- mycol = ['darkgreen', 'mediumaquamarine', 'silver', 'darkblue','steelblue']
- def cm2inch(value):
- return value/2.54
- def mt_examples(patient, save_path, fit_method, kmin, kmax,
- windowsize, windowstep, dt):
- recordings = mr_utilities.list_preictrecordings(patient)[0]
- fig, ax = plt.subplots(2, 1, figsize=(cm2inch(13.5), cm2inch(5)))
- binnings = ['right', 'left']
- quantity = 'm'
- t_onset = 600 # seizure onset at the end of each 10 min recording
- for irec, rec in enumerate(recordings):
- for ibin, binning in enumerate(binnings):
- mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_' \
- 'winstep{}_dt{}.pkl'.format(save_path, rec, binning,
- fit_method, kmin, kmax, windowsize, windowstep, dt)
- if not os.path.isfile(mt_file):
- print('{} not existing'.format(mt_file))
- continue
- # -----------------------------------------------------------------#
- # read m(t) data
- # -----------------------------------------------------------------#
- mt_frame = pd.read_pickle(mt_file)
- mt_frame = mr_utilities.exclude_inconsistencies(mt_frame)
- if mt_frame.empty:
- continue
- kmin = mt_frame.kmin.unique()[0]
- kmax = mt_frame.kmax.unique()[0]
- fit_method = mt_frame.fit_method.unique()[0]
- winsize = mt_frame.windowsize.unique()[0]
- t0_list = mt_frame.t0.tolist()
- qt = mt_frame[quantity].tolist()
- t0_list = [t0_ms * 0.001 for t0_ms in t0_list]
- half_win_sec = int((winsize / 2) / 1000 * dt)
- t_middle = [t0 + half_win_sec for t0 in t0_list]
- # -----------------------------------------------------------------#
- # Plotting
- # -----------------------------------------------------------------#
- ax[ibin].plot(t_middle, qt, '-', linewidth=1,
- label='Rec {}'.format(irec+1), color=mycol[irec])
- ax[ibin].axhline(y=1, color='black', linestyle='--', alpha=0.3)
- ax[ibin].set_ylabel(r"$\hat{m}$")
- ax[ibin].set_xlim(
- (t0_list[0] - 10, t0_list[-1] + 2 * half_win_sec + 10))
- binning_focus = mr_utilities.binning_focus(binning, int(patient))
- if binning_focus == 'focal':
- focus_label = 'ipsi'
- else:
- focus_label = 'contra'
- ax[ibin].set_title("{}".format(focus_label), loc='left')
- ax[0].legend(loc='center left', bbox_to_anchor=(1.04, -0.2),
- fontsize='small')
- ax[0].xaxis.grid(False)
- ax[1].xaxis.grid(False)
- ax[1].set_xlabel("time (s)")
- ax[0].annotate("", xy=(t_onset, 1), xytext=(t_onset, 1.1),
- arrowprops=dict(arrowstyle="->", color='black'))
- ax[1].annotate("", xy=(t_onset, 1), xytext=(t_onset, 1.04),
- arrowprops=dict(arrowstyle="->", color='black'))
- ax[0].annotate("seizure\n onset", xy=(t_onset, 1),
- xytext=(t_onset, 1.03))
- sns.despine()
- plt.subplots_adjust(bottom=None, right=None,
- wspace=0.3, hspace=1.3, left=None, top=None)
- plt.savefig('{}/Fig2_mt_patient{}.pdf'.format(save_path, patient),
- bbox_inches='tight')
- plt.close()
- def boxplot_mx(result_path, x, ksteps, fit_method, Engel=False):
- # -----------------------------------------------------------------#
- # Read and filter data
- # -----------------------------------------------------------------#
- mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
- result_path, x, ksteps[0], ksteps[1], fit_method)
- mx_allresults = pd.read_pickle(mx_results)
- mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults)
- x = mx_allresults['partno'].max()
- mx_allresults['logepsilon'] = mx_allresults.apply(lambda row: np.log(
- 1-row['m']), axis=1)
- mx_allresults['epsilon'] = mx_allresults.apply(lambda row: 1-row['m'],
- axis=1)
- measure = 'epsilon'
- patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
- if Engel:
- mx_allresults = mx_allresults[(mx_allresults['patient'].isin(patients1a))]
- # -----------------------------------------------------------------#
- # write estimates to lists
- # -----------------------------------------------------------------#
- focal_diff_list = []
- contralateral_diff_list = []
- all_last_contra = []
- all_first_contra = []
- all_last_focal = []
- all_first_focal = []
- for ir, recording in enumerate(mx_allresults['recording'].unique()):
- mx_recording = mx_allresults[mx_allresults['recording'] == recording]
- mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal']
- mx_contralateral = mx_recording[
- mx_recording['binning_focus'] == 'contra-lateral']
- if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[
- mx_focal['partno'] == x].empty:
- m_focal_first = mx_focal[mx_focal['partno'] == 1][
- measure].item()
- m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
- focal_diff_log = np.log(m_focal_last) - np.log(m_focal_first)
- focal_diff_list.append(focal_diff_log)
- all_first_focal.append(m_focal_first)
- all_last_focal.append(m_focal_last)
- if not mx_contralateral[mx_contralateral['partno'] == 1].empty \
- and not mx_contralateral[
- mx_contralateral['partno'] == x].empty:
- m_contralateral_first = \
- mx_contralateral[mx_contralateral['partno'] == 1][measure].item()
- m_contralateral_last = \
- mx_contralateral[mx_contralateral['partno'] == x][measure].item()
- contralateral_diff_log = np.log(m_contralateral_last) - \
- np.log(m_contralateral_first)
- contralateral_diff_list.append(contralateral_diff_log)
- all_first_contra.append(m_contralateral_first)
- all_last_contra.append(m_contralateral_last)
- # -----------------------------------------------------------------#
- # Wilcoxon paired test
- # -----------------------------------------------------------------#
- stat, p_contra_w = stats.wilcoxon(contralateral_diff_list)
- stat, p_focal_w = stats.wilcoxon(focal_diff_list)
- # -----------------------------------------------------------------#
- # Create dataframe for plotting
- # -----------------------------------------------------------------#
- quantity = measure
- data1 = {
- quantity: all_first_focal,
- 'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_first_focal))]
- }
- df = pd.DataFrame(data1)
- data2 = pd.DataFrame({
- quantity: all_last_focal,
- 'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
- })
- df = df.append(data2, ignore_index=True)
- data3 = pd.DataFrame({
- quantity: all_first_contra,
- 'rec_type': ['10 - 5 min' for _ in range(len(all_first_contra))],
- 'binning_focus': ['nSOZ' for _ in range(len(all_first_contra))]
- })
- df = df.append(data3, ignore_index=True)
- data4 = pd.DataFrame({
- quantity: all_last_contra,
- 'rec_type': ['5 - 0 min' for _ in range(len(all_last_contra))],
- 'binning_focus': ['nSOZ' for _ in range(len(all_last_contra))]
- })
- df = df.append(data4, ignore_index=True)
- # -----------------------------------------------------------------#
- # Plotting
- # -----------------------------------------------------------------#
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(5.5), cm2inch(4.5)))
- g = sns.boxplot(ax=ax, x='binning_focus', y=quantity, hue='rec_type',
- data=df, hue_order = ['10 - 5 min', '5 - 0 min'], order=[
- 'SOZ', 'nSOZ'], dodge=True, color='slategray', palette=[
- 'lightgray', 'slategray'], linewidth=1.1, width=0.6)
- g.set_yscale("log")
- plt.grid(axis='x', b=False)
- plt.rcParams['legend.title_fontsize'] = 9
- plt.legend(loc="upper left", ncol=1,
- fontsize='x-small', bbox_to_anchor=(0, 1.7), title='time to '
- 'seizure onset')
- n_f_pre = len(all_first_focal)
- n_f_inter = len(all_last_focal)
- n_c_pre = len(all_first_contra)
- n_c_inter = len(all_last_contra)
- ax.set_xticklabels(['ipsi\n({}, {})'.format(n_f_pre, n_f_inter),
- 'contra\n({}, {})'.format(n_c_pre, n_c_inter)])
- # -----------------------------------------------------------------#
- # quantity labels
- # -----------------------------------------------------------------#
- quantity_label = r'$\hat{m}$'
- ymax = 1e-3
- ymin = 1
- ax.set_ylim(ymin, ymax)
- T = [1, 1e-1, 1e-2, 1e-3]
- ax.set_yticks(T)
- ax.set_yticklabels([0, 0.9, 0.99, 0.999])
- # -----------------------------------------------------------------#
- # annotate p-values
- # -----------------------------------------------------------------#
- y_max = 2.4*ymax
- ax.annotate("", xy=(-0.2, y_max), xycoords='data',
- xytext=(0.2, y_max), textcoords='data',
- arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.2"))
- ax.text(-0.25, y_max*0.6, 'p={:.3f}'.format(p_focal_w), fontsize=8)
- ax.annotate("", xy=(0.8, y_max), xycoords='data',
- xytext=(1.2, y_max), textcoords='data',
- arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.2"))
- ax.text(0.75, y_max*0.6, 'p={:.3f}'.format(p_contra_w), fontsize=8)
- ax.set_ylim(ymin, ymax)
- plt.ylabel(quantity_label)
- plt.xlabel('')
- sns.despine()
- # -----------------------------------------------------------------#
- # Saving
- # -----------------------------------------------------------------#
- plotfile_mx = '{}Fig2_2parts_Engel{}.pdf'.format(result_path, Engel)
- plt.savefig(plotfile_mx, bbox_inches='tight')
- def pointplot_mx_regions_Engel(result_path, x, ksteps, fit_method):
- # -----------------------------------------------------------------#
- # Read and filter data
- # -----------------------------------------------------------------#
- mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
- result_path, x, ksteps[0], ksteps[1], fit_method)
- mx_allresults = pd.read_pickle(mx_results)
- mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults)
- x = mx_allresults['partno'].max()
- mx_allresults['epsilon'] = mx_allresults.apply(lambda row: 1-row['m'],
- axis=1)
- measure = 'epsilon'
- def give_region(row):
- return row['binning'][1:]
- mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row),
- axis=1)
- # -----------------------------------------------------------------#
- # write paired estimates to lists
- # -----------------------------------------------------------------#
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(4*1.2)))
- df = None
- patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
- mx_allresults = mx_allresults[(mx_allresults['patient'].isin(patients1a))]
- regions = ['H', 'A', 'PHC', 'EC']
- labels = []
- for i, reg in enumerate(regions):
- mx_region = mx_allresults[mx_allresults['region'] == reg]
- recordings = mx_region.recording.unique()
- all_first_focal = []
- all_last_focal = []
- all_first_contra = []
- all_last_contra = []
- recID = []
- for recording in recordings:
- mx_recording = mx_region[mx_region['recording'] == recording]
- mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal']
- mx_contralateral = mx_recording[
- mx_recording['binning_focus'] == 'contra-lateral']
- if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[
- mx_focal['partno'] == x].empty:
- m_focal_first = mx_focal[mx_focal['partno'] == 1][
- measure].item()
- m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
- all_first_focal.append(m_focal_first)
- all_last_focal.append(m_focal_last)
- recID.append(recording)
- if not mx_contralateral[mx_contralateral['partno'] == 1].empty \
- and not \
- mx_contralateral[
- mx_contralateral['partno'] == x].empty:
- m_contralateral_first = mx_contralateral[mx_contralateral[
- 'partno'] == 1][measure].item()
- m_contralateral_last = mx_contralateral[mx_contralateral[
- 'partno'] == x][measure].item()
- all_first_contra.append(m_contralateral_first)
- all_last_contra.append(m_contralateral_last)
- # -----------------------------------------------------------------#
- # Create dataframe for plotting
- # -----------------------------------------------------------------#
- quantity = measure
- if df is None:
- data1 = {
- 'region': reg,
- quantity: all_first_focal,
- 'recording': recID,
- 'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_first_focal))]
- }
- df = pd.DataFrame(data1)
- else:
- data1 = pd.DataFrame({
- 'region': reg,
- quantity: all_first_focal,
- 'recording': recID,
- 'rec_type': ['10 - 5 min' for _ in range(len(all_last_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
- })
- df = df.append(data1, ignore_index=True)
- data2 = pd.DataFrame({
- 'region': reg,
- quantity: all_last_focal,
- 'recording': recID,
- 'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
- })
- df = df.append(data2, ignore_index=True)
- n_f_first = len(all_first_focal)
- n_f_last = len(all_last_focal)
- labels.append('{}\n({}, {})'.format(reg, n_f_first, n_f_last))
- x0 = (-0.2+i) * np.ones(len(all_first_focal))
- x1 = (0.2+i) * np.ones(len(all_last_focal))
- plt.plot([x0, x1], [all_first_focal, all_last_focal],
- color='grey', linewidth=0.5, linestyle='--')
- # -----------------------------------------------------------------#
- # Plotting
- # -----------------------------------------------------------------#
- df = df[df['binning_focus'] == 'SOZ']
- g = sns.stripplot(ax=ax, x='region', y=quantity, hue='rec_type',
- data=df, hue_order=['10 - 5 min', '5 - 0 min'],
- dodge=True, order=regions, jitter=False,
- color='slategray', palette=['slategray', 'darkblue'],
- size=2.5, alpha=0.7)
- ax.set_yscale("log")
- ax.grid(axis='x', b=False)
- plt.rcParams['legend.title_fontsize'] = 9
- plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),
- fontsize='x-small', title='time to seizure onset')
- ax.set_xticklabels(labels)
- # -----------------------------------------------------------------#
- # quantity labels
- # -----------------------------------------------------------------#
- quantity_label = r'$\hat{m}$'
- ymax = 1e-3
- ymin = 1
- ax.set_ylim(ymin, ymax)
- T = [1, 1e-1, 1e-2, 1e-3]
- ax.set_yticks(T)
- ax.set_yticklabels([0, 0.9, 0.99, 0.999])
- # -----------------------------------------------------------------#
- # Annotate p-values of pairwise comparison
- # -----------------------------------------------------------------#
- ax.set_ylim(ymin, ymax)
- plt.ylabel(quantity_label)
- plt.xlabel('brain area')
- sns.despine()
- # -----------------------------------------------------------------#
- # Saving
- # -----------------------------------------------------------------#
- plotfile_mx = '{}Fig2_2parts_regions_Engel1a.pdf'.format(result_path)
- plt.savefig(plotfile_mx, bbox_inches='tight')
- def boxplot_mx_regions(result_path, x, ksteps, fit_method, Engel=False):
- # -----------------------------------------------------------------#
- # Read and filter data
- # -----------------------------------------------------------------#
- mx_results = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
- result_path, x, ksteps[0], ksteps[1], fit_method)
- mx_allresults = pd.read_pickle(mx_results)
- mx_allresults = mr_utilities.exclude_inconsistencies(mx_allresults)
- x = mx_allresults['partno'].max()
- mx_allresults['epsilon'] = mx_allresults.apply(lambda row:
- 1-row['m'], axis=1)
- measure = 'epsilon'
- def give_region(row):
- return row['binning'][1:]
- mx_allresults['region'] = mx_allresults.apply(lambda row: give_region(row),
- axis=1)
- if Engel:
- patients1a = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
- mx_allresults = mx_allresults[
- (mx_allresults['patient'].isin(patients1a))]
- # -----------------------------------------------------------------#
- # write estimates to lists
- # -----------------------------------------------------------------#
- df = None
- regions = ['H', 'A', 'PHC', 'EC']
- labels = []
- pvals = dict()
- for reg in regions:
- mx_region = mx_allresults[mx_allresults['region'] == reg]
- recordings = mx_region.recording.unique()
- focal_diff_list = []
- contralateral_diff_list = []
- all_first_focal = []
- all_last_focal = []
- all_first_contra = []
- all_last_contra = []
- recID = []
- for recording in recordings:
- mx_recording = mx_region[mx_region['recording'] == recording]
- mx_focal = mx_recording[mx_recording['binning_focus'] == 'focal']
- mx_contralateral = mx_recording[
- mx_recording['binning_focus'] == 'contra-lateral']
- if not mx_focal[mx_focal['partno'] == 1].empty and not mx_focal[
- mx_focal['partno'] == x].empty:
- m_focal_first = mx_focal[mx_focal['partno'] == 1][
- measure].item()
- m_focal_last = mx_focal[mx_focal['partno'] == x][measure].item()
- focal_diff_log = np.log(m_focal_last) - np.log(m_focal_first)
- all_first_focal.append(m_focal_first)
- all_last_focal.append(m_focal_last)
- focal_diff_list.append(focal_diff_log)
- recID.append(recording)
- if not mx_contralateral[mx_contralateral['partno'] == 1].empty \
- and not mx_contralateral[mx_contralateral['partno'] == x].empty:
- m_contralateral_first = mx_contralateral[mx_contralateral['partno'] == 1][measure].item()
- m_contralateral_last = mx_contralateral[mx_contralateral['partno'] == x][measure].item()
- contralateral_diff_log = np.log(m_contralateral_last) - \
- np.log(m_contralateral_first)
- all_first_contra.append(m_contralateral_first)
- all_last_contra.append(m_contralateral_last)
- contralateral_diff_list.append(contralateral_diff_log)
- # -----------------------------------------------------------------#
- # Wilcoxon paired test
- # -----------------------------------------------------------------#
- stat, p_contra_w = stats.wilcoxon(contralateral_diff_list)
- stat, p_focal_w = stats.wilcoxon(focal_diff_list)
- pvals[reg] = p_focal_w
- # -----------------------------------------------------------------#
- # Create dataframe for plotting
- # -----------------------------------------------------------------#
- quantity = measure
- if df is None:
- data1 = {
- 'region': reg,
- quantity: all_first_focal,
- 'recording': recID,
- 'rec_type': ['10 - 5 min' for _ in range(len(all_first_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_first_focal))]
- }
- df = pd.DataFrame(data1)
- else:
- data1 = pd.DataFrame({
- 'region': reg,
- quantity: all_first_focal,
- 'recording': recID,
- 'rec_type': ['10 - 5 min' for _ in range(len(all_last_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
- })
- df = df.append(data1, ignore_index=True)
- data2 = pd.DataFrame({
- 'region': reg,
- quantity: all_last_focal,
- 'recording': recID,
- 'rec_type': ['5 - 0 min' for _ in range(len(all_last_focal))],
- 'binning_focus': ['SOZ' for _ in range(len(all_last_focal))]
- })
- df = df.append(data2, ignore_index=True)
- data3 = pd.DataFrame({
- 'region': reg,
- quantity: all_first_contra,
- 'recording': ['-' for _ in range(len(all_first_contra))],
- 'rec_type': ['10 - 5 min' for _ in range(len(all_first_contra))],
- 'binning_focus': ['nSOZ' for _ in range(len(all_first_contra))]
- })
- df = df.append(data3, ignore_index=True)
- data4 = pd.DataFrame({
- 'region': reg,
- quantity: all_last_contra,
- 'recording': ['-' for _ in range(len(all_first_contra))],
- 'rec_type': ['5 - 0 min' for _ in range(len(all_last_contra))],
- 'binning_focus': ['nSOZ' for _ in range(len(all_last_contra))]
- })
- df = df.append(data4, ignore_index=True)
- n_f_first = len(all_first_focal)
- n_f_last = len(all_last_focal)
- labels.append('{}\n({}, {})'.format(reg, n_f_first, n_f_last))
- # -----------------------------------------------------------------#
- # Plotting
- # -----------------------------------------------------------------#
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(10*1.2), cm2inch(4*1.2)))
- df = df[df['binning_focus'] == 'SOZ']
- g = sns.boxplot(ax=ax, x='region', y=quantity, hue='rec_type',
- data=df, hue_order=['10 - 5 min', '5 - 0 min'],
- order=regions, dodge=True, color='slategray', palette=[
- 'lightgray','slategray'], linewidth=1.1, width=0.8)
- g.set_yscale("log")
- ax.grid(axis='x', b=False)
- plt.rcParams['legend.title_fontsize'] = 9
- plt.legend(loc='center left', bbox_to_anchor=(1, 0.5),
- fontsize='x-small', title='time to seizure onset')
- ax.set_xticklabels(labels)
- # -----------------------------------------------------------------#
- # quantity labels
- # -----------------------------------------------------------------#
- quantity_label = r'$\hat{m}$'
- ymax = 1e-3
- ymin = 1
- ax.set_ylim(ymin, ymax)
- T = [1, 1e-1, 1e-2, 1e-3]
- ax.set_yticks(T)
- ax.set_yticklabels([0, 0.9, 0.99, 0.999])
- # -----------------------------------------------------------------#
- # Annotate p-values of pairwise comparison
- # -----------------------------------------------------------------#
- y_max = 1.3 * ymax
- y_min = ymin
- for ir, region in enumerate(regions):
- if not np.isnan(pvals[region]):
- ax.annotate("", xy=(-0.2 + ir * 1.0, 2.0 * y_max), xycoords='data',
- xytext=(0.2 + ir * 1.0, 2.0 * y_max),
- textcoords='data',
- arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.3"))
- ax.text(-0.35 + ir * 1.0, 0.9*y_max * 1. - abs(y_max - y_min) *
- 0.0, 'p' + '={:.3f}'.format(pvals[region]), fontsize=8.5)
- ax.set_ylim(ymin, ymax)
- plt.ylabel(quantity_label)
- plt.xlabel('brain area')
- sns.despine()
- # -----------------------------------------------------------------#
- # Saving
- # -----------------------------------------------------------------#
- plotfile_mx = '{}Fig2_2parts_regions_Engel{}.pdf'.format(result_path,
- Engel)
- plt.savefig(plotfile_mx, bbox_inches='tight')
- def binning_label(patient, soz):
- focifile = '../Data/patients.txt'
- f = open(focifile, 'r')
- foci = [line.rstrip('\n').split('\t') for line in f.readlines()]
- focus = 'NA'
- for i, idf in enumerate(foci[1:]):
- if int(idf[0]) == int(patient):
- focus = idf[1]
- if soz == "SOZ":
- if focus == "L":
- binning = "left"
- elif focus == "R":
- binning = "right"
- else:
- raise Exception("No clear focus found.")
- elif soz == "nSOZ":
- if focus == "L":
- binning = "right"
- elif focus == "R":
- binning = "left"
- else:
- raise Exception("No clear focus found.")
- return binning
- def get_mt(patient, recording, soz, save_path, fit_method, kmin, kmax,
- windowsize, windowstep, dt):
- binning = binning_label(patient, soz)
- mt_file = '{}{}/mt_results_{}_{}_kmin{}_kmax{}_winsize{}_winstep{}_dt{}.pkl'.format(
- save_path, recording, binning, fit_method, kmin, kmax,
- windowsize, windowstep, dt)
- if not os.path.isfile(mt_file):
- raise Exception('{} not existing'.format(mt_file))
- # -----------------------------------------------------------------#
- # read m(t) data
- # -----------------------------------------------------------------#
- mt_frame = pd.read_pickle(mt_file)
- original_length = len(mt_frame.m.tolist())
- rejected_inconsistencies = ['H_nonzero', 'H_linear', 'not_converged']
- for incon in rejected_inconsistencies:
- if not mt_frame.empty:
- mt_frame = mt_frame[[incon not in mt_frame[
- 'inconsistencies'].tolist()[i] for i in range(
- len(mt_frame['inconsistencies'].tolist()))]]
- if len(mt_frame.m.tolist()) < 0.95 * original_length or len(
- mt_frame.m.tolist()) < 200: # in general, the recording should be
- # 260 windows long. if it is much shorter, it is unclear where to
- # divide it.
- raise Exception('Too many inconsistencies')
- mt_list = mt_frame.m.tolist()
- return mt_list
- def mad(data):
- return np.median(np.abs(data-np.median(data)))
- def variance_mt(save_path, data_path, fit_method, kmin, kmax,
- windowsize, windowstep, dt, logepsilon=True, Engel=False):
- patients = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
- 15, 16, 17, 18, 19, 20]
- if Engel:
- patients = [1, 6, 8, 9, 10, 11, 16, 17, 18, 19]
- laterality = ["SOZ", "nSOZ"]
- # -----------------------------------------------------------------#
- # Extract m(t) values an compute variances
- # -----------------------------------------------------------------#
- df = None
- for j, soz in enumerate(laterality):
- variance_first = []
- variance_last = []
- diff_list = []
- for ir, patient in enumerate(patients):
- recordings = mr_utilities.list_preictrecordings(patient)[0]
- for irec, rec in enumerate(recordings):
- try:
- m_list = get_mt(patient, rec, soz, data_path, fit_method,
- kmin, kmax, windowsize, windowstep, dt)
- m_first = m_list[:int(len(m_list)/2)]
- m_last = m_list[-int(len(m_list) / 2):]
- if logepsilon:
- m_first = [np.log(1 - m) for m in m_first]
- m_last = [np.log(1 - m) for m in m_last]
- variance_first.append(mad(m_first))
- variance_last.append(mad(m_last))
- diff_list.append(mad(m_first)-mad(m_last))
- except Exception:
- # If too many window estimates were excluded in the
- # function get_mt(), an exception will be risen.
- # We will simply continue with the remaining recordings.
- continue
- # -----------------------------------------------------------------#
- # create boxplot of variance
- # -----------------------------------------------------------------#
- stat, p = stats.wilcoxon(diff_list)
- print('p-value:'+str(p))
- quantity = 'variance'
- data1 = pd.DataFrame({
- quantity: variance_first,
- 'rec_type': ['10 - 5 min' for _ in range(len(variance_first))],
- 'binning_focus': [soz for _ in range(len(variance_first))],
- 'pval': [p for _ in range(len(variance_last))]
- })
- if df is None:
- df = data1
- else:
- df = df.append(data1, ignore_index=True)
- data2 = pd.DataFrame({
- quantity: variance_last,
- 'rec_type': ['5 - 0 min' for _ in range(len(variance_last))],
- 'binning_focus': [soz for _ in range(len(variance_last))],
- 'pval': [p for _ in range(len(variance_last))]
- })
- df = df.append(data2, ignore_index=True)
- # -----------------------------------------------------------------#
- # Plotting
- # -----------------------------------------------------------------#
- fig, ax = plt.subplots(1, 1, figsize=(cm2inch(5.5), cm2inch(4.5)))
- g = sns.boxplot(ax=ax, x='binning_focus', y='variance', hue='rec_type',
- data=df, hue_order=['10 - 5 min', '5 - 0 min'], order=[
- 'SOZ', 'nSOZ'], dodge=True, color='slategray', palette=[
- 'lightgray', 'slategray'], linewidth=1.1, width=0.6)
- g.set_yscale("log")
- plt.grid(axis='x', b=False)
- plt.rcParams['legend.title_fontsize'] = 9
- plt.legend(loc="upper left", ncol=1, fontsize='x-small',
- bbox_to_anchor=(0, 1.7), title='time to seizure onset')
- n_f_pre = len(df[(df['binning_focus'] == 'SOZ') &
- (df['rec_type'] == '5 - 0 min')].index)
- n_f_inter = len(df[(df['binning_focus'] == 'SOZ') &
- (df['rec_type'] == '10 - 5 min')].index)
- n_c_pre = len(df[(df['binning_focus'] == 'nSOZ') &
- (df['rec_type'] == '5 - 0 min')].index)
- n_c_inter = len(df[(df['binning_focus'] == 'nSOZ') &
- (df['rec_type'] == '10 - 5 min')].index)
- ax.set_xticklabels(['ipsi\n({}, {})'.format(n_f_pre, n_f_inter),
- 'contra\n({}, {})'.format(n_c_pre, n_c_inter)])
- # -----------------------------------------------------------------#
- # annotate p-values
- # -----------------------------------------------------------------#
- y_max = 2.4 * np.max(variance_first)
- y_min = 0.4 * np.min(variance_first)
- ax.annotate("", xy=(-0.2, y_max), xycoords='data',
- xytext=(0.2, y_max), textcoords='data',
- arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.2"))
- ax.text(-0.2, y_max + abs(y_max - y_min) * 0.1,
- 'p={:.3f}'.format(df[df['binning_focus'] == 'SOZ'].pval.tolist()[
- 0]), fontsize=8)
- ax.annotate("", xy=(0.8, y_max), xycoords='data', xytext=(1.2, y_max),
- textcoords='data', arrowprops=dict(arrowstyle="-", ec='black',
- connectionstyle="bar,fraction=0.2"))
- ax.text(0.8, y_max + abs(y_max - y_min) * 0.1,
- 'p={:.3f}'.format(df[df['binning_focus'] == 'nSOZ'].pval.tolist(
- )[0]), fontsize=8)
- plt.xlabel('')
- sns.despine()
- plt.ylabel(r'Variability in $\hat{m}$')
- plt.yscale('log')
- plt.savefig('{}Fig2_variance_mt.pdf'.format(save_path), bbox_inches='tight')
- plt.show()
- def timeresolved_example():
- result_path = '../Results/preictal/singleunits/'
- ksteps = (1, 400)
- kmin = ksteps[0]
- kmax = ksteps[1]
- windowsize = 20000
- windowstep = 500
- dt = 4
- patients = [13]
- for patient in patients:
- mt_examples(patient, result_path, 'offset', kmin, kmax,
- windowsize, windowstep, dt)
- def boxplot_2parts_hemispheres():
- x = 2
- ksteps = (1, 400)
- fit_method = 'offset'
- result_path = '../Results/preictal/singleunits/'
- boxplot_mx(result_path, x, ksteps, fit_method, Engel=False)
- boxplot_mx(result_path, x, ksteps, fit_method, Engel=True)
- def boxplot_2parts_regions():
- x = 2
- ksteps = (1, 400)
- fit_method = 'offset'
- result_path = '../Results/preictal/singleunits_regions/'
- boxplot_mx_regions(result_path, x, ksteps, fit_method, Engel=False)
- pointplot_mx_regions_Engel(result_path, x, ksteps, fit_method)
- def visualize_mt_variance():
- save_path = '../Results/preictal/singleunits/'
- data_path = '../Results/preictal/singleunits/'
- variance_mt(save_path, data_path, fit_method='offset', kmin=1, kmax=400,
- windowsize=20000, windowstep=500, dt=4,
- Engel=False, logepsilon=True)
- def main():
- timeresolved_example()
- boxplot_2parts_hemispheres()
- boxplot_2parts_regions()
- visualize_mt_variance()
- if __name__ == '__main__':
- main()
|