diff --git a/.DS_Store b/.DS_Store index 2c9f02d..211e1bc 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/berlin_constants.py b/berlin_constants.py index dfcb4a5..4e141c3 100644 --- a/berlin_constants.py +++ b/berlin_constants.py @@ -91,161 +91,7 @@ ], } layout = BIDSLayout(data_path) -files_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - session=['EcogLfpMedOn01', 'EcogLfpMedOn02', 'EcogLfpMedOn03'], - return_type="filename", - ) -files1_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='001', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files1_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='003', - session=['EcogLfpMedOn01', 'EcogLfpMedOn02', 'EcogLfpMedOn03'], - return_type="filename", - ) -files3_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='003', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files3_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='003', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files4_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='004', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files4_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='004', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files5_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='005', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files5_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='005', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files6_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='006', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files6_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='006', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files7_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='007', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files7_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='007', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files8_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='008', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files8_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='008', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files9_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='009', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files9_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='009', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files10_off = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='010', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) -files10_on = layout.get( - extension="vhdr", - task="Rest", - acquisition="StimOff", - subject='010', - session=['EcogLfpMedOff01','EcogLfpMedOff02', 'EcogLfpMedOff03'], - return_type="filename", - ) \ No newline at end of file +files = layout.get( + extension="vhdr", task="Rest", acquisition="StimOff", return_type="filename", +) + diff --git a/berlin_pipeline.py b/berlin_pipeline.py index e9fde9e..2e83ce8 100644 --- a/berlin_pipeline.py +++ b/berlin_pipeline.py @@ -1,10 +1,12 @@ +from typing import Union import runpy +from turtle import done import pandas as pd from bids import BIDSLayout import mne_bids - +import matplotlib.pyplot as plt from src import pipeline @@ -17,41 +19,48 @@ def main(): path_bids = project_constants["PATH_BIDS"] m1_ids = project_constants["M1_IDS"] new_ch_names_map = project_constants["NEW_CH_NAMES_MAP"] - files_off = project_constants["files_off"] - files_on = project_constants["files_on"] - files1_off = project_constants["files1_off"] - files1_on = project_constants["files1_on"] - files3_off = project_constants["files3_off"] - files3_on = project_constants["files3_on"] - files4_off = project_constants["files4_off"] - files4_on = project_constants["files4_on"] - files5_off = project_constants["files5_off"] - files5_on = project_constants["files5_on"] - files6_off = project_constants["files6_off"] - files6_on = project_constants["files6_on"] - files7_off = project_constants["files7_off"] - files7_on = project_constants["files7_on"] - files8_off = project_constants["files8_off"] - files8_on = project_constants["files8_on"] - files9_off = project_constants["files9_off"] - files9_on = project_constants["files9_on"] - files10_off = project_constants["files10_off"] - files10_on = project_constants["files10_on"] + files = project_constants["files"] + remove_subjects: Union[str, None] = ["002", "010", "011", "012"] + if remove_subjects: + for remove_subject in remove_subjects: + files = [file for file in files if remove_subject not in file] + + files1_off = [f for f in files if "001" in f and "MedOff" in f] + files1_on = [f for f in files if "001" in f and "MedOn" in f] + files3_off = [f for f in files if "003" in f and "MedOff" in f] + files3_on = [f for f in files if "003" in f and "MedOn" in f] + files4_off = [f for f in files if "004" in f and "MedOff" in f] + files4_on = [f for f in files if "004" in f and "MedOn" in f] + files5_off = [f for f in files if "005" in f and "MedOff" in f] + files5_on = [f for f in files if "005" in f and "MedOn" in f] + files6_off = [f for f in files if "006" in f and "MedOff" in f] + files6_on = [f for f in files if "006" in f and "MedOn" in f] + files7_off = [f for f in files if "007" in f and "MedOff" in f] + files7_on = [f for f in files if "007" in f and "MedOn" in f] + files8_off = [f for f in files if "008" in f and "MedOff" in f] + files8_on = [f for f in files if "008" in f and "MedOn" in f] + files9_off = [f for f in files if "009" in f and "MedOff" in f] + files9_on = [f for f in files if "009" in f and "MedOn" in f] + files10_off = [f for f in files if "010" in f and "MedOff" in f] + files10_on = [f for f in files if "010" in f and "MedOn" in f] # Define variables burst_char_pd_all = [] M1_burst_dynamics_all = [] npow_list_all = [] - # Process Files # - for path_run in files3_off: + burst_char_pd_all_on = [] + M1_burst_dynamics_all_on = [] + npow_list_all_on = [] + + # Process runs in one subject # + # def bursts_single_subject(runs): + for path_run in files: entities = mne_bids.get_entities_from_fname(path_run) sub = entities["subject"] - ( - burst_char_pd, - M1_burst_dynamics, - npow, - ) = pipeline.bursts_single_subject( + med = "On" if "MedOn" in entities["session"] else "Off" + run = entities["run"] + (burst_char_pd, M1_burst_dynamics, npow,) = pipeline.bursts_single_run( path_run=path_run, path_bids=path_bids, sub=sub, @@ -59,21 +68,17 @@ def main(): new_ch_names=new_ch_names_map[sub], ) burst_char_pd["Subject"] = sub + burst_char_pd["Medication"] = med + burst_char_pd["Run"] = run M1_burst_dynamics["Subject"] = sub npow["Subject"] = sub burst_char_pd_all.append(burst_char_pd) M1_burst_dynamics_all.append(M1_burst_dynamics) npow_list_all.append(npow) -# burst_char_pd_df = pd.DataFrame(burst_char_pd_all) -# M1_burst_dynamics_df = pd.DataFrame(M1_burst_dynamics_all) -# npow_df = pd.DataFrame(npow_list_all) - - - return burst_char_pd_all, M1_burst_dynamics_all, npow_list_all if __name__ == "__main__": burst_char_pd_all, M1_burst_dynamics_all, npow_list_all = main() - print("done") +print(done) diff --git a/src/IO.py b/src/IO.py index fa0bb30..04dd840 100644 --- a/src/IO.py +++ b/src/IO.py @@ -1,6 +1,7 @@ import mne_bids import numpy as np + def read_BIDS_data(PATH_RUN, BIDS_PATH): """Given a run path and bids data path, read the respective data Parameters @@ -17,20 +18,20 @@ def read_BIDS_data(PATH_RUN, BIDS_PATH): entities = mne_bids.get_entities_from_fname(PATH_RUN) bids_path = mne_bids.BIDSPath( - subject=entities["subject"], - session=entities["session"], - task=entities["task"], - run=entities["run"], - acquisition=entities["acquisition"], - datatype="ieeg", - root=BIDS_PATH, - ) + subject=entities["subject"], + session=entities["session"], + task=entities["task"], + run=entities["run"], + acquisition=entities["acquisition"], + datatype="ieeg", + root=BIDS_PATH, + ) raw_arr = mne_bids.read_raw_bids(bids_path) return ( raw_arr, raw_arr.get_data(), int(np.ceil(raw_arr.info["sfreq"])), - int(raw_arr.info["line_freq"]),) - + int(raw_arr.info["line_freq"]), + ) diff --git a/src/__pycache__/pipeline.cpython-39.pyc b/src/__pycache__/pipeline.cpython-39.pyc index fbc0720..4587d94 100644 Binary files a/src/__pycache__/pipeline.cpython-39.pyc and b/src/__pycache__/pipeline.cpython-39.pyc differ diff --git a/src/__pycache__/preprocessing.cpython-39.pyc b/src/__pycache__/preprocessing.cpython-39.pyc index 492c5cf..b39a2c3 100644 Binary files a/src/__pycache__/preprocessing.cpython-39.pyc and b/src/__pycache__/preprocessing.cpython-39.pyc differ diff --git a/src/burst_calc.py b/src/burst_calc.py index 6573132..ff8d8a3 100644 --- a/src/burst_calc.py +++ b/src/burst_calc.py @@ -9,42 +9,48 @@ def Time_Frequency_Estimation(signal): - freqs = np.arange(1,101) - power = mne.decoding.TimeFrequency(freqs, sfreq= 1600 , method='morlet', n_cycles=10,decim=8, output='power' ) + freqs = np.arange(1, 101) + power = mne.decoding.TimeFrequency( + freqs, sfreq=1600, method="morlet", n_cycles=10, decim=8, output="power" + ) run_TF = power.transform(signal) - return (run_TF) + return run_TF + def beta_bands(run_TF): - ''' + """ Beta bands of the ecog channels: low beta(13-20Hz), high beta (20-35Hz), full beta (13-35Hz) - ''' + """ l_low_beta = [] l_high_beta = [] l_full_beta = [] for ch_idx in range(run_TF.shape[0]): - l_low_beta.append(run_TF[ch_idx, LOW_BETA[0]:LOW_BETA[1],:]) - l_high_beta.append(run_TF[ch_idx, HIGH_BETA[0]:HIGH_BETA[1],:]) - l_full_beta.append(run_TF[ch_idx, FULL_BETA[0]:FULL_BETA[1],:]) + l_low_beta.append(run_TF[ch_idx, LOW_BETA[0] : LOW_BETA[1], :]) + l_high_beta.append(run_TF[ch_idx, HIGH_BETA[0] : HIGH_BETA[1], :]) + l_full_beta.append(run_TF[ch_idx, FULL_BETA[0] : FULL_BETA[1], :]) return l_low_beta, l_high_beta, l_full_beta + def avg_power(l_beta): return [np.mean(l_ch, axis=0) for l_ch in l_beta] + def z_score(l_beta): return [stats.zscore(l_ch, axis=0) for l_ch in l_beta] + def percentile(l_beta, percentile): return [np.percentile(l_ch, q=percentile) for l_ch in l_beta] -def get_burst_length(beta_averp_norm,beta_thr, sfreq=200): +def get_burst_length(beta_averp_norm, beta_thr, sfreq=200): """ Analysing the duration of beta burst """ bursts = np.zeros((beta_averp_norm.shape[0] + 1), dtype=bool) bursts[1:] = beta_averp_norm >= beta_thr - deriv = np.diff(bursts) + deriv = np.diff(bursts) isburst = False burst_length = [] burst_start = 0 @@ -60,21 +66,22 @@ def get_burst_length(beta_averp_norm,beta_thr, sfreq=200): isburst = True if isburst: burst_length.append(index + 1 - burst_start) - burst_length = np.array(burst_length)/sfreq - + burst_length = np.array(burst_length) / sfreq + return burst_length + def exclude_short_bursts(burst_length): - ''' + """ exclude bursts shorter than 100ms - ''' - return [i for i in burst_length if i >= 0.1] + """ + return [i for i in burst_length if i >= 0.1] def get_burst_amplitude(beta_amplitude, beta_thr): - amplitude = [] + amplitude = [] no_bursts = True - cont = False + cont = False for val in beta_amplitude: if val >= beta_thr: no_bursts = False @@ -92,9 +99,3 @@ def get_burst_amplitude(beta_amplitude, beta_thr): return mean_amplitude - - - - - - diff --git a/src/pipeline.py b/src/pipeline.py index 0d25da0..ecc0bd5 100644 --- a/src/pipeline.py +++ b/src/pipeline.py @@ -9,7 +9,7 @@ from src import burst_calc, IO, postprocessing, preprocessing -def bursts_single_subject( +def bursts_single_run( path_run: Union[str, pathlib.Path], path_bids: Union[str, pathlib.Path], sub: str, @@ -21,13 +21,9 @@ def bursts_single_subject( raw_ecog = preprocessing.pick_ecog(raw) if sub == "001": - raw_ecog_bi = preprocessing.bipolar_reference_s1( - raw, raw_ecog, new_ch_names - ) + raw_ecog_bi = preprocessing.bipolar_reference_s1(raw, raw_ecog, new_ch_names) else: - raw_ecog_bi = preprocessing.bipolar_reference( - raw, raw_ecog, new_ch_names - ) + raw_ecog_bi = preprocessing.bipolar_reference(raw, raw_ecog, new_ch_names) NUM_CH = len(raw_ecog_bi.get_channel_types()) @@ -51,9 +47,7 @@ def bursts_single_subject( l_beta_avg_norm = [burst_calc.z_score(l) for l in l_beta_avg] # 75th percentile of the power - l_beta_thr = [ - burst_calc.percentile(l, percentile=75) for l in l_beta_avg_norm - ] + l_beta_thr = [burst_calc.percentile(l, percentile=75) for l in l_beta_avg_norm] low = 0 high = 1 @@ -83,16 +77,13 @@ def bursts_single_subject( burst_duration_cl_m1 = burst_duration_cl[m1] # Mean burst duration mean_burst_duration = [ - np.nanmean(burst_duration_cl[ch_idx], axis=0) - for ch_idx in range(NUM_CH) + np.nanmean(burst_duration_cl[ch_idx], axis=0) for ch_idx in range(NUM_CH) ] mean_dur_m1 = mean_burst_duration[m1] # Histogram of burst duration histogram_duration = [ - np.histogram( - burst_duration_cl[ch_idx], density=False, bins=20, range=(0, 2) - )[0] + np.histogram(burst_duration_cl[ch_idx], density=False, bins=20, range=(0, 2))[0] for ch_idx in range(NUM_CH) ] hist_dur_m1 = histogram_duration[m1] @@ -128,9 +119,7 @@ def bursts_single_subject( ) # M1 Beta Burst Dynamics - M1_burst_dynamics = postprocessing.dataframe_burst_dynamics( - normhist_dur_m1 - ) + M1_burst_dynamics = postprocessing.dataframe_burst_dynamics(normhist_dur_m1) # normalized beta power npow = postprocessing.dataframe_npow(psd_M1) diff --git a/src/plot_utils.py b/src/plot_utils.py index 08ad65a..7dd6399 100644 --- a/src/plot_utils.py +++ b/src/plot_utils.py @@ -6,7 +6,7 @@ ALPHA_BOX = 0.4 -def histplot_burst_length(burst_length, bins=13, range=(0, 1.3)): +def histplot_burst_length(burst_length, bins=9, range=(0, 0.9)): """ plot in histogram burst lengths""" sns.set(style="white", font_scale=1) plt.hist(burst_length, bins=bins, range=range) diff --git a/src/preprocessing.py b/src/preprocessing.py index 21f9ac0..3c6fe72 100644 --- a/src/preprocessing.py +++ b/src/preprocessing.py @@ -2,68 +2,83 @@ import mne from scipy import stats + def pick_ecog(raw): - ''' + """ pick ECoG channels - ''' + """ raw_ecog = [] raw_ecog.append(raw.pick_types(ecog=True).ch_names) raw_ecog_red = raw_ecog[0][0:6] return raw_ecog_red -def bipolar_reference (raw, raw_ecog, new_ch_names): - + +def bipolar_reference(raw, raw_ecog, new_ch_names): + anode = raw_ecog[0:5] cathode = raw_ecog[1:6] - raw_ecog_bi = mne.set_bipolar_reference(raw.load_data(), anode=anode, - cathode=cathode, ch_name=new_ch_names ) - return (raw_ecog_bi) + raw_ecog_bi = mne.set_bipolar_reference( + raw.load_data(), anode=anode, cathode=cathode, ch_name=new_ch_names + ) + return raw_ecog_bi + + +def bipolar_reference_s1(raw, raw_ecog, new_ch_names): -def bipolar_reference_s1 (raw, raw_ecog, new_ch_names): - anode = raw_ecog[0:4] cathode = raw_ecog[1:5] - raw_ecog_bi = mne.set_bipolar_reference(raw.load_data(), anode=anode, cathode=cathode, ch_name=new_ch_names ) + raw_ecog_bi = mne.set_bipolar_reference( + raw.load_data(), anode=anode, cathode=cathode, ch_name=new_ch_names + ) raw_ecog_bi = raw_ecog_bi.load_data().add_channels( + [raw.pick_channels(["ECOG_L_1_2_SMC_AT"])], force_update_info=True + ) + channel_order = [ + "ECOG_L_1_2_SMC_AT", + "ECOG_L_2_3_SMC_AT", + "ECOG_L_3_4_SMC_AT", + "ECOG_L_4_5_SMC_AT", + "ECOG_L_5_6_SMC_AT", + ] + raw_ecog_bi = raw_ecog_bi.reorder_channels(channel_order) + return raw_ecog_bi + - [raw.pick_channels(["ECOG_L_1_2_SMC_AT"])], - force_update_info=True -) - channel_order = ['ECOG_L_1_2_SMC_AT','ECOG_L_2_3_SMC_AT', - 'ECOG_L_3_4_SMC_AT', - 'ECOG_L_4_5_SMC_AT', - 'ECOG_L_5_6_SMC_AT'] - raw_ecog_bi = raw_ecog_bi.reorder_channels(channel_order) - return (raw_ecog_bi) - -def bipolar_reference_s10_on (raw, raw_ecog, new_ch_names): - +def bipolar_reference_s10_on(raw, raw_ecog, new_ch_names): anode = raw_ecog[0:4] cathode = raw_ecog[1:5] - raw_ecog_bi = mne.set_bipolar_reference(raw.load_data(), anode=anode, - cathode=cathode, ch_name=new_ch_names ) - return (raw_ecog_bi) + raw_ecog_bi = mne.set_bipolar_reference( + raw.load_data(), anode=anode, cathode=cathode, ch_name=new_ch_names + ) + return raw_ecog_bi -def filtering (raw_ecog_bi): - ''' + +def filtering(raw_ecog_bi): + """ Filtering (Highpass 3 Hz, Notchfilter 50,251,50Hz, Lowpass 250Hz) - ''' + """ raw_ecog_hi = raw_ecog_bi.copy().filter(3, None,) raw_ecog_hi_lo = raw_ecog_hi.copy().filter(None, 250) - raw_ecog_filt = raw_ecog_hi_lo.copy().notch_filter(np.arange(50,251,50), filter_length="auto", phase='zero') - return(raw_ecog_filt) + raw_ecog_filt = raw_ecog_hi_lo.copy().notch_filter( + np.arange(50, 251, 50), filter_length="auto", phase="zero" + ) + return raw_ecog_filt + def downsample(raw_ecog_filt): - ''' + """ Downsample Data to 1600Hz for faster processing - ''' + """ raw_ecog_dow = raw_ecog_filt.copy().resample(1600) - return (raw_ecog_dow) + return raw_ecog_dow + -def get_data (raw_ecog_dow): +def get_data(raw_ecog_dow): signal = raw_ecog_dow.copy().get_data() - return (signal) + return signal + def z_score_signal(signal): stand_signal = stats.zscore(signal, axis=1) - return (stand_signal) \ No newline at end of file + return stand_signal +