# Import librairies
import random
import neo.rawio
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
import seaborn as sb
import os
import neo
import mne
import pandas as pd
import scipy.stats as stats
import matplotlib
[docs]
def ensure_dir(path: str) -> None:
"""
Create the directory if it does not exist
Parameters
---------------------------
path: string
Path-like where you want to create folder
"""
# Check if the path exists
isExist = os.path.exists(path)
# If not create it
if not isExist:
os.makedirs(path)
[docs]
def load_raw_data(path,dtype,length,pattern2exclude, time='Random',analog=False,*args):
"""
Load data from the raw files (e.g., .ncs for Neuralynx)
This function supports multiple file formats:
- .ncs for Neuralynx acquisition system
- .nsX for Blaclrock acquisition system
- .med for Darkhorseneuro acquisition system
- .edf
This function is based on MNE [1] and NEO [2] framework
Returns two objects :
1) A raw-object from MNE (see https://mne.tools/stable/generated/mne.io.Raw.html#mne.io.Raw for documentation)
2) An array containing raw values of the signal you want to analyze
Parameters
---------------------------
path: string
A string containing the path where your data are stored
dtype: string
For Neuralynx data : 'ncs'
For Blackrock data : 'nsX'
For Dark Horse Neuro data : 'med'
edf is also supported : 'edf'
length: int or 'all'
If you specify length with a int (e.g., 5) it will select randomly 5 minute of signal
If you specify 'all' it will take the entire recording (can be slower, depends on your computational resources)
pattern2exclude: string
Pattern that enables us to not take into account the macrocontacts
(e.g., in Toulouse all macro contacts have the suffix _sub so the argument will be '_sub.ncs')
time: string or tuple (Default='Random')
By default take randomly portion of signal corresponding to length
It can be a tuple (2,7) for example to take signal from 2 minutes to 7 minutes of the recording
analog: string (optional) (Default: False)
Reference channel used in blackrock recording system if it exists
Returns
---------------------------
raw_data: raw data in FIF format (see MNE)
This object contains all informations about your data (channels name, sampling rate etc.)
data_sig: np array
Array containing raw values of the signal you want to load
Will have a shape of nCh x nSamples
References
----------
[1] Gramfort A, Luessi M, Larson E, Engemann DA, Strohmeier D, Brodbeck C, Parkkonen L, Hämäläinen MS.
MNE software for processing MEG and EEG data.
Neuroimage. 2014 Feb 1;86:446-60.
doi: 10.1016/j.neuroimage.2013.10.027.
Epub 2013 Oct 24. PMID: 24161808; PMCID: PMC3930851.
[2] Garcia S., Guarino D., Jaillet F., Jennings T.R., Pröpper R., Rautenberg P.L.,
Rodgers C., Sobolev A.,Wachtler T., Yger P. and Davison A.P. (2014)
Neo: an object model for handling electrophysiology data in multiple formats.
Frontiers in Neuroinformatics 8:10: doi:10.3389/fninf.2014.00010
"""
# Load 'ncs' from Neuralynx
if dtype=='ncs':
# Exclude files containing 'sub' because correspond to macro-electrodes - TO adapt with a pattern in arguments of the fucntions
raw_data = mne.io.read_raw_neuralynx(path,exclude_fname_patterns=list(['*'+pattern2exclude]))
# Load 'ns5' from Blackrock
elif dtype=='nsX':
raw_data = mne.io.read_raw_nsx(path)
if analog:
# Exclude the reference channel
raw_data.drop_channels(analog)
# Load edf data files
elif dtype=='edf':
raw_data = mne.io.read_raw_edf(path)
# If there is a trigger channels in the edf we want to exclude it
try:
raw_data.drop_channels('trigger')
except:
print('No trigger channels to drop')
# Pas sur qu'on doive garder ça....
elif dtype=='dat':
# Here not sure everyone is using int16 but for us it is ok
data_type = np.int16
# In order to read the dat we absolutely need the number of channels so if not specified by user return a message
if len(args)==0:
print("Please indicate the number of channels in your dat file to load it")
# Here we load data
elif len(args)==1:
size = os.path.getsize(path)
size = int(size/np.dtype(data_type).itemsize)
raw_data = np.memmap(path, mode='r', dtype=data_type, order='F', shape=(args[0], int(size/args[0])))
# If the user specify too much arguments we ask him to only add the number the channels and nothing more
else:
print("You specify too much arguments please only indicate the number of channels")
elif dtype == 'med':
# Here we create a raw object from the med folder using this:
# https://mne.tools/stable/auto_examples/io/read_neo_format.html
# Read data from the .med folder
reader = neo.io.MedIO(path)
# Get the first block - proxy (do not load in memory)
block = reader.read(lazy=True)[0]
# Get a proxy data from first segment (you have to get only one segment)
segment = block.segments[0]
# Get the signal from all channels
signals_proxy = segment.analogsignals[0]
else:
print("File format not supported, you can load .ncs, .nsX, .med, .edf for now")
print("Contact us to add new file format")
print("dtype must be either ncs, nsX, med, edf")
# If you want to load all recording
if length=='all':
# If want to load all file then load all med recording
if dtype == 'med':
signals = signals_proxy.load()
data = signals.rescale("V").magnitude.T
sfreq = signals.sampling_rate.magnitude
# Get the name of the channels
ch_names = [f"{reader.header['signal_channels'][idx][0]}" for idx in range(signals.shape[1])]
# Attribute a type to channels (here eeg)
ch_types = ["eeg"] * len(ch_names) # if not specified, type 'misc' is assumed
# Create a raw object in MNE
info = mne.create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sfreq)
raw_data = mne.io.RawArray(data, info)
# Specify tmin=0 so it starts at the beginning of the recording
raw_interest = raw_data.crop(tmin=0)
# Load data into memory
data_sig = raw_interest.get_data()
# For all other file format easier to load data
else:
raw_interest = raw_data.crop(tmin=0)
data_sig = raw_interest.get_data()
# Else will randomly choosed a portion of length specified in input
else:
if dtype == 'med':
import random
# Determine the index max for the length you want
# e.g. If you want 5 minutes of signal the last sample can not be less than 5 minutes before the end of your recording
sampling_rate = int(signals_proxy.sampling_rate)
last_sample = int(signals_proxy.duration*sampling_rate)
idx_max_random = int(last_sample - (sampling_rate*60*length)) # (sampling_rate*60*length) correspond to the length of your subselection
# Here randomly select the starting index of the 5 minutes
if time == 'Random':
idx_time = int(random.random()*idx_max_random)
signal_crop = signals_proxy.time_slice(t_start=idx_time/sampling_rate,t_stop=(idx_time/sampling_rate)+length*60)
else:
# Crop signal from the portion of interest
signal_crop = signals_proxy.time_slice(t_start=time[0]*60,t_stop=(time[1]*60)) # transform time in seconds
data = signal_crop.rescale("V").magnitude.T
sfreq = signal_crop.sampling_rate.magnitude
# Get the name of the channels
ch_names = [f"{reader.header['signal_channels'][idx][0]}" for idx in range(signal_crop.shape[1])]
# Attribute a type to channels (here eeg)
ch_types = ["eeg"] * len(ch_names) # if not specified, type 'misc' is assumed
# Create a raw object in MNE
info = mne.create_info(ch_names=ch_names, ch_types=ch_types, sfreq=sfreq)
raw_data = mne.io.RawArray(data, info)
# Load data in memory
data_sig = raw_data.get_data()
else:
if time == 'Random':
# Select randomly the first sample
idx_time = random_time(raw_data,raw_data.info['sfreq'],length)
sampling_rate = raw_data.info['sfreq']
# Crop the signal between 1st sample and last sample (last sample - 1st sample = length)
raw_interest = raw_data.crop(tmin=idx_time/sampling_rate,tmax=(idx_time/sampling_rate)+length*60)
else:
raw_interest = raw_data.crop(tmin=time[0]*60,tmax=time[1]*60)
# Load data in memory
data_sig = raw_interest.get_data()
return raw_data,data_sig
[docs]
def splitCharNum(string) :
"""
Separate a string containing characters and numbers in two objects
(One with characters and one with numbers)
Parameters
---------------------------
string: string
String containing characters and numbers (e.g., 'da1')
Returns
---------------------------
char: string
string containing only the characters from the input (e.g., 'da')
number: string
string containing only the numbers from the input (e.g., '1')
"""
import re
char,number=re.findall(r'[A-Za-z-_\'\s]+|\d+', string)
return char,number
[docs]
def get_unique_unsorted(array):
"""
Maintains the order of appearance in the original array
Does not sort by ascending order (figures) or alphabetic order (characters)
Parameters
---------------------------
array: array of string
Array of string containing the name of your channels (e.g. ['da','db'])
Returns
---------------------------
uniqueUnsorted: string
array of string containing
"""
unique, uniqueInds=np.unique(array, return_index=True)
uniqueUnsorted=array[np.sort(uniqueInds)]
return uniqueUnsorted
[docs]
def reorder_data(Regions,data):
'''
When there is 3 tetrodes on the same shaft order will be x1, x10, x11, x12, x2, x3 ... x9 on the raw file
We want to re-organize it so the order is x1, x2, x3 ... x9, x10, x11, x12
Will return the array of regions name and data in the right order
Parameters
---------------------------
Regions: array of string
Array of strings containing all the name of the regions implanted with your channels
data: array
Matrix containing raw values with shape nChannels x nSamples
Returns
---------------------------
Regions_ok: array of string
Array of strings containing the name in the right order (e.g., 'x1', 'x2', 'x3' ... 'x9', 'x10', 'x11', 'x12')
data_ok: array
Matrix with shape nChannels x nSamples but re-organize to match the order of channel labels
'''
# Extract the label of the electrode
chChar=np.array([splitCharNum(reg)[0] for reg in Regions])
# Extract the number of the channels recording
chNumber=np.array([splitCharNum(reg)[1] for reg in Regions])
# Get the index to put them in the right order
regInds=np.array([chChar==char for char in get_unique_unsorted(chChar)])
argsortChRegs=np.concatenate([np.where(regInds[regi])[0][np.argsort(np.array(chNumber[regInds[regi]], int))]\
for regi in range(regInds.shape[0])])
# Re-order the names and lfps in right order
Regions_ok = Regions[argsortChRegs]
data_ok = data[argsortChRegs]
return Regions_ok,data_ok
[docs]
def find_nearest(array, value):
"""
Find the index where there is the nearest values in an array from the one we want
Parameters
---------------------------
array: np.array
array containing data (e.g., [2, 7, 12])
value: int
The value we want to find the closest (e.g. 11)
Returns
---------------------------
idx: int
Will return the index in array closest to value (e.g., 2)
"""
# Make sure it is an array
array = np.asarray(array)
# Get the index where closest from value
idx = (np.abs(array - value)).argmin()
return idx
[docs]
def random_time(data,sampling_rate,length):
"""
Select randomly length minutes of signal from your entire recoding
Parameters
---------------------------
data: MNE-raw object
The structure of a RAW mne object containing metadata
sampling rate: int
The sampling rate of your recording (e.g., 32768)
length: int
Length in minute of the signal you want to analyze
Returns
---------------------------
start_idx: int
Correspond to the first sample of the signal, the beginning of the recording you want
to analyze
"""
# Determine the index max for the length you want
# e.g. If you want 5 minutes of signal the last sample can not be less than 5 minutes before the end of your recording
last_sample = data.last_samp.T
idx_max_random = int(last_sample - (sampling_rate*60*length)) # (sr*60*lenght) correspond to the length of your subselection
# Here randomly select the starting index of the 5 minutes
start_idx = int(random.random()*idx_max_random)
return start_idx
[docs]
def p_welch(data,chRegions,sub,sess,sr,sr_down,fr_low,fr_high,saveFolder,probe_type='Dixi'):
"""
Compute the power spectrum of your signal in order to identify frequencies present in your signal
Here we use the welch's method to compute the power spectrum.
Parameters
---------------------------
data: array
Matrice containing data of all channels containing nSamples
chRegions: Numpy array
Array containing the labels of all the electrodes in the recording
sub: string
ID of the patient you are analyzing
sess: string
Name of the session you are analyzing
sr: int
Sampling rate of your recording (e.g., 32768)
sr_down: int
Sampling rate downsampled, so the one you want (e.g., 8192 Hz)
fr_low: int
The lowest frequency from which we compute the power spectrum
fr_high: int
The highest frequency from which we compute the power spectrum
saveFolder: string
Path where you want to save the output figure
probe_type: string (Default=Dixi)
String containing the model of your electrodes, it can be either 'Dixi' or 'Ad-tech'
Returns
---------------------------
pxx_log: array
Array containing power spectrum values for the channel of interest
f_plot: array
Array containing the values of frequency associate with each power spectrum value (will be useful for the plot)
"""
# Downsample based on sr_down
if sr_down != 1:
data_ds = sig.decimate(data,sr_down)
# Store the down sampled data into data
data = data_ds
# Initialize the list where to store the power values
pxx_log = []
# Loop over all channels in your recording
for iCh in range(data.shape[0]):
# Compute welch method to estimate PSD
f, pxx = sig.welch(data[iCh], fs=sr, nperseg=4096)
# Get the index from fr_low to fr_high Hz to plot the results
idx_fr_lim = find_nearest(f,fr_high)
idx_debut = find_nearest(f,fr_low)
# Store values of the power spectrum
pxx_log.append(10*np.log10(pxx[idx_debut:idx_fr_lim]))
f_plot = f[idx_debut:idx_fr_lim]
# Do the plot part
plot_all_chan(f_plot,data.shape[0],chRegions,pxx_log,sub,sess,probe_type,saveFolder)
return pxx_log
[docs]
def plot_all_chan(f,nCh,chRegs,psd,bsnm,session,probe_type,saveFolder):
"""
This function will plot the power spectrum of all micro-channels on the same plot
Parameters
---------------------------
f: array
Array containing frequency values (it is the output of p_welch)
nCh: int
Number of channels in your recording
chRegs: array of strings
Array containing name of your channels
psd: array
Array of power spectrum values for each value of f, it is the output of p_welch
session: string
Name of the session you analyze, specify at the beginning of the jupyter notebook
probe_type: string
String containing the model of micro-electrodes you have in your dataset (for now: Dixi or Ad-tech only)
saveFolder: string
String containing the path of the folder where you want to save figure
Returns
---------------------------
Matplotlib plot containing power spectrum of each channels and saved in the folder specified
"""
# Make sure to create the directory
ensure_dir(saveFolder)
# List of colors, each tetrode will be in the same color but with different transparency
colors = ['black','red','orangered','saddlebrown','gold','olive','chartreuse','turquoise','darkslategray','dodgerblue','midnightblue','slateblue','darkviolet','violet','magenta','crimson']
# Get different alpha levels for each wire
if probe_type == 'Dixi':
transparency = [0.4,0.6,0.8,1]
elif probe_type == 'Ad-tech':
transparency = [0.4,0.5,0.6,0.7,0.8,0.9,1]
# Initialize to zero the current tetrode (so it takes the first one)
iGroup = 0
# Initialize the lists where we will store values
min_pwr = []
max_pwr = []
# Plot the PSD of the micro-wire with according alpha transparency
fig, ax1 = plt.subplots(1,1,layout='constrained')
# Here we loop on all channel and plot the PSD corresponding to each channel
for i in range(0,nCh):
# Determine wich micro-wire of the tetrode it is (1st, 2nd, 3rd, 4th)
if probe_type == 'Dixi':
modulo_ch = i%4
elif probe_type == 'Ad-tech':
modulo_ch = i%8
ax1.plot(f,psd[i],color=colors[iGroup],alpha=transparency[modulo_ch])
# Keep min and max to automatically adjust the limit of the plot
min_pwr.append(min(psd[i]))
max_pwr.append(max(psd[i]))
# When we did the last micro-wire of the tetrode iTetrode increment to go to the next tetrode
if probe_type == 'Dixi':
if modulo_ch==3:
iGroup=iGroup+1
elif probe_type == 'Ad-tech':
if modulo_ch == 7:
iGroup = iGroup + 1
# Legend and save plot
ax1.set_ylabel('10 * log10(Power)')
ax1.set_xlabel('Frequency (Hz)')
ax1.set_title('Power Spectrum (Welch ''s method) - '+bsnm+' - '+session)
ax1.legend(loc='center left',labels=chRegs,bbox_to_anchor=(1.01,0.5),fontsize=4)
# Save the figure in the right folder
plt.rcParams["svg.fonttype"] = 'none'
plt.savefig(saveFolder + 'PSD_All_Channels_'+bsnm+'_'+session+'.svg', dpi=300)
# Display the figure in the notebook
plt.show()
# Close it
plt.close()
[docs]
def plot_tetrode(metrics,chRegions):
'''
Function used to plot the metrics value of each tetrodes grouped by colors
'''
# List of colors, each tetrode will be in the same color but with different transparency
colors = ['black','red','orangered','saddlebrown','gold','olive','chartreuse','turquoise','darkslategray','dodgerblue','midnightblue','slateblue','darkviolet','violet','magenta','crimson']
# Loop over all tetrodes in your recording
for iTetrode in range(int(len(chRegions)/4)):
# Extract the metrics of the tetrode analyzed
metrics_tetrode = [metrics[iTetrode*4],metrics[(iTetrode*4)+1],metrics[(iTetrode*4)+2],metrics[(iTetrode*4)+3]]
# Plot each tetrode with the same color
plt.plot(chRegions[iTetrode*4:iTetrode*4+4],metrics_tetrode,color=colors[iTetrode])
[docs]
def plot_noise(data,sr,chRegions,path,save=1,limit='auto',fr_low=300,fr_high=3000):
'''
This function will plot the raw signal filtered between 300 and 300 Hz for 1 second randomly choosed in the 5 minutes window
So we can have an idea of the level of noise during our recording
Parameters
---------------------------
data: array
Matrix with your data with shape nChannels x nSamples
sr: int
sampling rate of the signal
chRegions: string
Array of string containing name of each channels
path: string
String containing the path where you want to save the figures
save: int, default = 1
If you don't want to save figure put save = 0
limit: string
Default value is 'auto' so matplotlib handes the limit but otherwise it is the min and max of your data
Could be useful to better see the level of noise in your recording
fr_low: int, default = 300
The lowest frequency for your band pass filter
fr_high: int, default = 3000
The highest frequency for your band pass filter
Returns
---------------------------
Matplotlib plot saved in the specific folder if you want to
'''
# Check whether the specified path exists or not
ensure_dir(path)
# Select 1s of data
idx_max_random = int(data.shape[1] - (sr)) # Select the index so the beginning is at least 1s before the end of the recording
# Here randomly select the starting index of the second
idx_debut = int(random.random()*idx_max_random)
# Design our butterworth
try:
sos = sig.butter(3,[fr_low,fr_high],'bandpass',fs=sr,output='sos')
except:
sos = sig.butter(3,[fr_low,(sr/2)-1],'bandpass',fs=sr,output='sos')
# Initialize the list where to store data filtered between 300 and 3000 Hz
data_filtered = list()
# Loop over all the channels
for iCh in range(data.shape[0]):
# Filter the data and store them in the list
data_filtered.append(sig.sosfilt(sos,data[iCh]))
# Get the max value
limit_min = np.min(data_filtered)
limit_max = np.max(data_filtered)
# Plot results
for iCh in range(data.shape[0]):
# Plot the data
fig, ax1 = plt.subplots(1,1,layout='constrained',figsize=(12,5))
ax1.plot(np.linspace(0,1,sr),data_filtered[iCh][idx_debut:idx_debut+sr]*1000000,color='cornflowerblue')
ax1.set_ylabel('µV',color='grey',fontsize=15)
ax1.tick_params(axis='y',colors='grey',labelsize=12)
ax1.set_xlabel('Time (s)',fontsize=15)
ax1.tick_params(axis='x',labelsize=12)
ax1.set_title('Noise level - channel : '+ chRegions[iCh] + ' (n° : '+str(iCh)+')',fontsize=18)
# Remove axis lines on top and right part of the box
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
if limit != 'auto':
ax1.set_ylim((limit_min,limit_max))
if save==1:
plt.savefig(path+'Noise_level_channel_'+chRegions[iCh]+'.jpg')
plt.close()
[docs]
def plot_raw(data,sr,num_channel,path,save=1):
'''
In entry take the 5 minutes of signal that were randomly choosed before
This function will plot the raw signal for 1 second randomly choosed in the 5 minutes window
Parameters
---------------------------
data: array
data from one particular channel
sr: int
sampling rate of the signal
num_channel: int
Number of the channel in your recording
path: string
Path where you want to store your figure in output
save: Boolean
Default = 1 else put save = 0
Returns
---------------------------
Matplotlib plot with raw signal
'''
# Check whether the specified path exists or not
ensure_dir(path)
# Select 1s of data
idx_max_random = int(len(data) - (sr)) # Select the index so the beginning is at least 1s before the end of the recording
# Here randomly select the starting index of the second
idx_debut = int(random.random()*idx_max_random)
# Plot the data
plt.figure(figsize=(12,5))
plt.plot(np.linspace(0,1,sr),data[idx_debut:idx_debut+sr])
plt.ylabel('µV')
plt.xlabel('Time (s)')
plt.title('Raw signal - channel : ' + str(num_channel))
if save==1:
plt.savefig(path+'Raw_Signal_channel'+str(num_channel)+'.jpg')
plt.close()
[docs]
def tblprep(path,electrodes,sub,sess) :
'''
Here we prep the excel file where all metrics values will be stored
You can either create a file for each session or patient by changing the name of path
But you could also keep the same file and then store values to create a database.
Parameters
---------------------------
path: string
Path where you want to store the excel file containing the metrics or an existing file to store values
electrodes: array of string
Array containing the name of your electrodes (automatically extracted from your recording)
sub: string
Name of the subject analyzed
sess: string
Name of the session analyzed
Outputs
---------------------------
It creates an excel file where metrics will be stored
'''
#open the preexisting excel file
import os.path as path_os
# It is an existing excel load it
if path_os.exists(path):
tbl = pd.read_excel(path, header=0)
else:
# if empty, creat columns with named sub, session, and with electrode names
tbl = pd.DataFrame()
tbl['sub'] = 'DefaultValue'
tbl['run'] = 'DefaultValue'
tbl['electrodes'] = 'DefaultValue'
tbl['RMS_filter'] = 'DefaultValue'
tbl['variance'] = 'DefaultValue'
tbl['variance_norm'] = 'DefaultValue'
tbl['tetrode_cor'] = 'DefaultValue'
tbl['deviation'] = 'DefaultValue'
tbl['kurtosis'] ='DefaultValue'
tbl['region'] = 'DefaultValue'
tbl['SNR'] = 'DefaultValue'
tbl['Artefact'] = 'DefaultValue'
tbl['Hurst'] = 'DefaultValue'
#find the first empty line
if len(tbl) == 0:
i=0
else :
i=0
while i < len(tbl) and pd.notna(tbl.iloc[i, 0]) and tbl.iloc[i, 0] != "":
i = i + 1
#write sub, session and electrode names if they not already exist
if len(tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess)]) == 0:
for j in range(i,i+len(electrodes)):
x = j - i
tbl = pd.concat([tbl, pd.DataFrame([{'sub': sub, 'run': sess,'electrodes': electrodes[x]}])], ignore_index=True)
#save and replace the old file
tbl.to_excel(path, index=False)
[docs]
def rms_signal_filtered (data,path,chRegions,sub,sess,fr_low,fr_high,sr,saveFolder,save=1) :
'''
Here we take the 5 minutes of signal that were randomly choosed before
This function will calculate the RMS (root mean square) for 1 second randomly choosed in the 5 minutes window
Parameters
---------------------------
data: numpy.array
2-D Matrice containing all your data with shape nChannels x nSamples
path: string
Path where is the excel file you created earlier
chRegions: array of string
Name of the electrodes in your recordings
sub: string
Name of the subject analyzed
sess: string
Name of the session analyzed
fr_low: int
Low frequency to filter the data (e.g., 300)
fr_high: int
High frequency to filter the data (e.g., 3000)
sr: int
Sampling rate of the signal of your recording system
saveFolder: string
Path where you want to store the outputes figure
Outputs
---------------------------
Figure with the values of RMS of your signal filtered
'''
# Band-pass filter between two frequencies
# Design our filter
sos = sig.butter(3,[fr_low,fr_high],'bandpass',fs=sr,output='sos')
# Initialize the list
rms = list()
# Loop over the channels in the recording
for iChannels in range(int(data.shape[0])):
# Filter the data
data_filtered = sig.sosfilt(sos,data[iChannels])
# Calculate the Root Mean Square (RMS) on the whole data
rms.append(np.sqrt(np.mean(data_filtered**2)))
#open the table
tbl = pd.read_excel(path)
#write in the table
#write the correlation in the table
#creat a datframe with electrode names and rms values
rms_df = pd.DataFrame()
rms_df['electrodes'] = chRegions
rms_df['RMS_filter'] = rms
#write the correlation in the table
for i in range (0, len(rms_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == rms_df.iloc[i, 0]), 'RMS_filter'] = rms_df.iloc[i, 1]
#save the table
tbl.to_excel(path, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5),layout='constrained')
plot_tetrode(rms,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.title('RMS (filtered signal) for each channel',fontsize=18)
plt.xticks(rotation = 45)
plt.ylabel('RMS (300 - 3000 Hz) [Volts]',fontsize=15)
if save==1:
plt.savefig(saveFolder+'RMS_Filter_AllMicro.jpg')
return rms
[docs]
def plot_rms_filter(pathtbl,savingpath,sub,sess,save=1):
'''
This function is used to plot the Root-mean square value of each channels
Parameters
---------------------------
'''
tbl = pd.read_excel(pathtbl)
rmstbl = tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess), ['sub', 'run', 'electrodes', 'RMS_filter']]
plt.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5))
plt.title('RMS (filtered signal) for each channel',fontsize=18)
plt.xticks(rotation = 45)
plt.ylabel('RMS (300 - 3000 Hz) [Volts]',fontsize=15)
#sb.lineplot(x='electrodes', y='RMS_filter', data=rmstbl)
RMS_filter = rmstbl['RMS_filter']
chRegions = rmstbl['electrodes']
plot_tetrode(RMS_filter,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
if save==1:
plt.savefig(savingpath+'RMS_Filter_AllMicro.jpg')
[docs]
def correlation_coefficient(data,chRegions,path,pathtbl,sub,sess,probe_type,save=1):
'''
This function compute a correlation coefficient between one micro-wire and all the other micro-wires
from the same tetrode
see Tuyisenge, V., Trebaul, L., Bhattacharjee, M., Chanteloup-Forêt, B., Saubat-Guigui, C., Mîndruţă, I., Rheims, S.,
Maillard, L., Kahane, P., Taussig, D., & David, O. (2018).
Automatic bad channel detection in intracranial electroencephalographic recordings using ensemble machine learning.
Clinical Neurophysiology, 129(3), 548‑554. https://doi.org/10.1016/j.clinph.2017.12.013
Parameters
---------------------------
data: array
Matrix with your data with shape nChannels x nSamples
chRegions: Array
Vector with the name of the micro-wires
path: string
Path where you want to save the figure obtained
pathtbl: string
Path where you store the excel file containing all values for the patient you analyze
sub: string
String like 'sub-XX' with XX being the number of the subject in your database
sess: string
String containing the name of the session you are analyzing (specified at the beginning of the jupyter notebook)
probe_type: string
String containing the model of micro-electrodes you have in your dataset (for now: Dixi or Ad-tech only)
save: int (default=1)
To save the figure, put save = 0 if you don't want to save the figure
Returns
---------------------------
Matplotlib plot with correlation coefficient for each channel
Values of the correlation coefficient saved in an excel fil to get a report of the analyzes
'''
import os
# Check whether the specified path exists or not
ensure_dir(path)
mean_corr = list()
# Loop over each tetrode
if probe_type == 'Dixi':
for iTetrode in range(int(data.shape[0]/4)):
# Select data for the tetrode of interest
data_tetrode = data[iTetrode*4:iTetrode*4+4]
# For each channel in the tetrode of interest
for i in range(data_tetrode.shape[0]):
corr = np.zeros(data_tetrode.shape[0])
for j in range(data_tetrode.shape[0]):
if i!=j:
# Compute the pearson correlation between channel i and j
corr[j] = stats.pearsonr(data[i+iTetrode*4],data[j+iTetrode*4]).statistic
# Get the mean correlation with all neighbouring channels
corr[corr==0] = np.nan
mean_corr.append(np.nanmean(corr))
elif probe_type == 'Ad-tech':
for iGroup in range(int(data.shape[0]/8)):
# Select data for the tetrode of interest
data_group = data[iGroup*8:iGroup*8+8]
# For each channel in the tetrode of interest
for i in range(data_group.shape[0]):
corr = np.zeros(data_group.shape[0])
for j in range(data_group.shape[0]):
if i!=j:
# Compute the pearson correlation between channel i and j
corr[j] = stats.pearsonr(data[i+iGroup*8],data[j+iGroup*8]).statistic
# Get the mean correlation with all neighbouring channels
corr[corr==0] = np.nan
mean_corr.append(np.nanmean(corr))
#open the excel table
tbl = pd.read_excel(pathtbl)
#creat a datframe with electrode names and correlation
cor_df = pd.DataFrame()
cor_df['electrodes'] = chRegions
cor_df['tetrode_cor'] = mean_corr
#write the correlation in the table
for i in range (0, len(cor_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == cor_df.iloc[i, 0]), 'tetrode_cor'] = cor_df.iloc[i, 1]
#save table
tbl.to_excel(pathtbl, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5),layout='constrained')
plot_tetrode(mean_corr,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.title('Correlation coefficient with neighbouring channels',fontsize=18)
plt.xticks(rotation = 45)
plt.ylabel('Pearson (r-value)',fontsize=15)
if save==1:
plt.savefig(path+'Correlation_Coefficient_AllMicroChannels.jpg')
return mean_corr
[docs]
def variance_normalized(data,chRegions,path,pathtbl,sub,sess,probe_type,save=1):
'''
Compute the variance of each channel normalized by the mean variance of all neighbouring channels on the same tetrode
If variance is high then it means that the channel is not recording the same activity than the ones around, so maybe it
is broken or there is a problem with the plugging on your recording system
Parameters
---------------------------
data: array
Matrix with your data with shape nChannels x nSamples
chRegions: array of strings
Vector with the name of your channels
path: string
Path where you want to save output figure
pathtbl: string
Path where the excel file of this patient is stored
sub: string
String containing the number or the ID of the patient in your database
sess: string
String containing the name of the session you are analyzing
probe_type: string
String containing the model of micro-electrodes you have in your dataset (for now: Dixi or Ad-tech only)
save: int (Default=1)
To save figure or not (put = 0 to not save)
Returns
---------------------------
Matplotlib plot with values of variance normalized
Values are also stored in the excel file
see Tuyisenge, V., Trebaul, L., Bhattacharjee, M., Chanteloup-Forêt, B., Saubat-Guigui, C., Mîndruţă, I., Rheims, S.,
Maillard, L., Kahane, P., Taussig, D., & David, O. (2018).
Automatic bad channel detection in intracranial electroencephalographic recordings using ensemble machine learning.
Clinical Neurophysiology, 129(3), 548‑554. https://doi.org/10.1016/j.clinph.2017.12.013
'''
# Check whether the specified path exists or not
ensure_dir(path)
variance = list()
if probe_type == 'Dixi':
# Loop over all tetrodes
for iTetrode in range(int(data.shape[0]/4)):
# Select data of interest
data_tetrode = data[iTetrode*4:iTetrode*4+4]
# Loop over all channels in the tetrode of interest
for i in range(data_tetrode.shape[0]):
# Get the variance of channel i
var_i = np.var(data_tetrode[i])
var_j = np.zeros(data_tetrode.shape[0])
for j in range(data_tetrode.shape[0]):
if i!=j:
# Get the variance of neighbouring channels
var_j[j] = np.var(data_tetrode[j])
# Get the mean variance of neighbouring electrodes
var_j[var_j==0] = np.nan
median_var_j = np.nanmedian(var_j)
# Normalized variance of channels i by mean variance of neighbouring channels
variance.append(var_i/median_var_j)
elif probe_type == 'Ad-tech':
# Loop over all tetrodes
for iGroup in range(int(data.shape[0]/8)):
# Select data of interest
data_group = data[iGroup*8:iGroup*8+8]
# Loop over all channels in the tetrode of interest
for i in range(data_group.shape[0]):
# Get the variance of channel i
var_i = np.var(data_group[i])
var_j = np.zeros(data_group.shape[0])
for j in range(data_group.shape[0]):
if i!=j:
# Get the variance of neighbouring channels
var_j[j] = np.var(data_group[j])
# Get the mean variance of neighbouring electrodes
var_j[var_j==0] = np.nan
mean_var_j = np.nanmean(var_j)
# Normalized variance of channels i by mean variance of neighbouring channels
variance.append(var_i/mean_var_j)
#open the excel table
tbl = pd.read_excel(pathtbl)
#creat a datframe with electrode names and
var_df = pd.DataFrame()
var_df['electrodes'] = chRegions
var_df['var'] = variance
#write the correlation in the table
for i in range (0, len(var_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == var_df.iloc[i, 0]), 'variance_norm'] = var_df.iloc[i, 1]
#save table
tbl.to_excel(pathtbl, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5),layout='constrained')
plot_tetrode(variance,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.title('Variance normalized by neighbouring mico-channels',fontsize=18)
plt.ylabel('Variance normalized',fontsize=15)
plt.xticks(rotation = 45)
if save==1:
plt.savefig(path+'Variance_Normalized_AllChannels.jpg')
[docs]
def deviation(data,chRegions,path,pathtbl,sub,sess,probe_type, save=1):
'''
Compute the deviation (i.e. electrical drift)
Parameters
---------------------------
data: array
Matrix with your data with shape nChannels x nSamples
chRegions: array of strings
Vector with the name of your channels
path: string
Path where you want to save the output's figure
paththbl : string
Path where to store values in excel
sub: string
IDs of the patient
sess: string
Name of the session
probe_type: string
String containing the model of micro-electrodes you have in your dataset (for now: Dixi or Ad-tech only)
Returns
---------------------------
Matplotlib plot and values in excel file
see Tuyisenge, V., Trebaul, L., Bhattacharjee, M., Chanteloup-Forêt, B., Saubat-Guigui, C., Mîndruţă, I., Rheims, S.,
Maillard, L., Kahane, P., Taussig, D., & David, O. (2018).
Automatic bad channel detection in intracranial electroencephalographic recordings using ensemble machine learning.
Clinical Neurophysiology, 129(3), 548‑554. https://doi.org/10.1016/j.clinph.2017.12.013
'''
# Check whether the specified path exists or not
ensure_dir(path)
deviation = list()
if probe_type == 'Dixi':
# Loop over all tetrodes
for iTetrode in range(int(data.shape[0]/4)):
# Select data of interest
data_tetrode = data[iTetrode*4:iTetrode*4+4]
# Loop over all channels in the tetrode of interest
for i in range(data_tetrode.shape[0]):
# Get the mean amplitude of channel i
mean_i = np.mean(data_tetrode[i])
mean_j = np.zeros(data_tetrode.shape[0])
for j in range(data_tetrode.shape[0]):
if i!=j:
# Get the mean amplitude of neighbouring channels
mean_j[j] = np.mean(data_tetrode[j])
# Get the mean of neighbouring electrodes' amplitudes
mean_j[mean_j==0] = np.nan
mean_neighbours = np.nanmean(mean_j)
# Get the deviation by substracting the mean of neighbours to channel i
deviation.append(mean_i - mean_neighbours)
elif probe_type == 'Ad-tech':
# Loop over all tetrodes
for iGroup in range(int(data.shape[0]/8)):
# Select data of interest
data_group = data[iGroup*8:iGroup*8+8]
# Loop over all channels in the tetrode of interest
for i in range(data_group.shape[0]):
# Get the mean amplitude of channel i
mean_i = np.mean(data_group[i])
mean_j = np.zeros(data_group.shape[0])
for j in range(data_group.shape[0]):
if i!=j:
# Get the mean amplitude of neighbouring channels
mean_j[j] = np.mean(data_group[j])
# Get the mean of neighbouring electrodes' amplitudes
mean_j[mean_j==0] = np.nan
mean_neighbours = np.nanmean(mean_j)
# Get the deviation by substracting the mean of neighbours to channel i
deviation.append(mean_i - mean_neighbours)
# Z-score transformation
deviation = stats.zscore(deviation)
#open the excel table
tbl = pd.read_excel(pathtbl)
#creat a datframe with electrode names and
dev_df = pd.DataFrame()
dev_df['electrodes'] = chRegions
dev_df['devi'] = deviation
#write the correlation in the table
for i in range (0, len(dev_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == dev_df.iloc[i, 0]), 'deviation'] = dev_df.iloc[i, 1]
#save table
tbl.to_excel(pathtbl, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5),layout='constrained')
plot_tetrode(deviation,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_title('Deviation - Electrical drift',fontsize=18)
ax.set_ylabel('Deviation',fontsize=15)
plt.xticks(rotation = 45)
if save==1:
plt.savefig(path+'Deviation_Zscore_AllChannels.jpg')
[docs]
def variance(data,chRegions,path,pathtbl,sub,sess,save=1):
'''
Compute the variance of each channel to see if there is a lot of artefacts or not
Parameters
---------------------------
data: array
Matrix with your data with shape nChannels x nSamples
chRegions: array of strings
Vector with the name of your channels
path: string
Path where you want to save output figure
paththbl: string
Path to the excel file where to store informations
sub: string
IDs of the patient
sess: string
Name of the session
save: int (Default=1)
To save figure
Returns
---------------------------
Matplotlib plot and value store in excel file
'''
# Check whether the specified path exists or not
ensure_dir(path)
variance = list()
# Loop over all tetrodes
for iChannels in range(int(data.shape[0])):
# Select data of interest
data_channel = data[iChannels]
# Get the variance of channel i
var_i = np.var(data_channel)
# Normalized variance of channels i by mean variance of neighbouring channels
variance.append(var_i)
#open the excel table
tbl = pd.read_excel(pathtbl)
#creat a datframe with electrode names and
var_df = pd.DataFrame()
var_df['electrodes'] = chRegions
var_df['var'] = variance
#write the correlation in the table
for i in range (0, len(var_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == var_df.iloc[i, 0]), 'variance'] = var_df.iloc[i, 1]
#save table
tbl.to_excel(pathtbl, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5),layout='constrained')
plot_tetrode(variance,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_title('Variance - Artifact',fontsize=18)
ax.set_ylabel('Variance',fontsize=15)
plt.xticks(rotation = 45)
if save==1:
plt.savefig(path+'Variance_AllChannels.jpg')
[docs]
def signaltonoise(a,chRegions,path,pathtbl,sub,sess, save=1, axis=1, ddof=0):
"""
The signal-to-noise ratio of the input data.
Returns the signal-to-noise ratio of `a`, here defined as the mean
divided by the standard deviation.
Parameters
---------------------------
a: array_like
An array_like object containing the sample data.
axis: int or None, optional
If axis is equal to None, the array is first ravel'd. If axis is an
integer, this is the axis over which to operate. Default is 0.
ddof: int, optional
Degrees of freedom correction for standard deviation. Default is 0.
Returns
---------------------------
s2n: ndarray
The mean to standard deviation ratio(s) along `axis`, or 0 where the
standard deviation is 0.
"""
# Check whether the specified path exists or not
ensure_dir(path)
a = np.asanyarray(a)
m = a.mean(axis)
sd = a.std(axis=axis, ddof=ddof)
signal2noise = np.where(sd == 0, 0, m/sd)
#open the excel table
tbl = pd.read_excel(pathtbl)
#creat a datframe with electrode names and
var_df = pd.DataFrame()
var_df['electrodes'] = chRegions
var_df['SNR'] = signal2noise
#write the correlation in the table
for i in range (0, len(var_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == var_df.iloc[i, 0]), 'SNR'] = var_df.iloc[i, 1]
#save table
tbl.to_excel(pathtbl, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5), layout='constrained')
plot_tetrode(signal2noise,chRegions)
plt.title('SNR')
plt.xticks(rotation = 45)
if save==1:
plt.savefig(path+'SNR_AllChannels.jpg')
return signal2noise
[docs]
def kurtosis(data,chRegions,path,pathtbl,sub,sess,save=1):
"""
Kurtosis: An electrical activity may appear in one of the channels and be absent in the remaining ones.
Such events can be detected by computing the kurtosis in all channels. Given that the kurtosis
indicates the presence of outliers in datasets, the highest value reveals which channel shows a
particular event (Mognon et al., 2011) (from Tuyisenge et al., 2018)
Parameters
---------------------------
data: array
Matrix containing your data with shape nChannels x nSamples
chRegions: array of strings
Vector containing name of your channels
path: string
String of the path where you want to store figures
pathtbl: string
Path where the excel file is stored
sub: string
IDs of the patient you are analyzing
sess: string
Name of the session
save: Boolean, default=1
=1 if you want to save plot, put = 0 otherwise
Results
---------------------------
Matplotlib plot and value store in the excel file
"""
kurtosis_channel = []
# Check whether the specified path exists or not
ensure_dir(path)
for iCh in range(data.shape[0]):
a = data[iCh]
# Compute kurtosis on each channel
kurtosis_channel.append(stats.kurtosis(a,axis=0,fisher=True))
#open the excel table
tbl = pd.read_excel(pathtbl)
#creat a datframe with electrode names and
kur_df = pd.DataFrame()
kur_df['electrodes'] = chRegions
kur_df['kurt'] = kurtosis_channel
#write the correlation in the table
for i in range (0, len(kur_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == kur_df.iloc[i, 0]), 'kurtosis'] = kur_df.iloc[i, 1]
#save table
tbl.to_excel(pathtbl, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5),layout='constrained')
plot_tetrode(kurtosis_channel,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_title('Kurtosis - Outliers in dataset',fontsize=18)
ax.set_ylabel('Kurtosis',fontsize=15)
plt.xticks(rotation = 45)
if save==1:
plt.savefig(path+'Kurtosis_AllChannels.jpg')
[docs]
def hurst_component(data,chRegions,path,pathtbl,sub,sess,save=1):
"""
Compute the Hurst component
You can see Tuyisenge et al. (2018) for a detail of the algorithm
Typically EEG data have values around 0.7
Parameters
---------------------------
data: array
Matrix containing your data with shape nChannels x nSamples
chRegions: array of strings
Vectors containing the name of your channels
path: string
Path where you want to store the plots
pathtbl: string
Path where is store your excel file saving results
sub: String
String "sub-XX" where XX is the number of the subject analyzed
sess: String
Name of the session you are analyzing
Returns
---------------------------
Matplotlib plot
"""
# Check whether the specified path exists or not
ensure_dir(path)
# Initialize list to store results
hurst = list()
# Loop over each channel
for iCh in range(data.shape[0]):
# Step 1. Compute mean amplitude
mean_amplitude = np.mean(data[iCh])
# Step 2. Create a mean centered channel
mean_center_channel = data[iCh] - mean_amplitude
# Step 3. Compute the cumulative channel deviation
chan_deviation = 0
chan_deviation = [chan_deviation + mean_center_channel[t] for t in range(len(mean_center_channel))]
# Step 4. Compute the channel amplitude range
amplitude_range = np.max(chan_deviation) - np.min(chan_deviation)
# Step 5. Compute the standard deviation
std_channel = np.std(data[iCh])
# Step 6. Compute Hurst exponent
hurst.append(np.sqrt(np.log(amplitude_range/std_channel)))
#open the excel table
tbl = pd.read_excel(pathtbl)
#creat a datframe with electrode names and
hurst_df = pd.DataFrame()
hurst_df['electrodes'] = chRegions
hurst_df['Hurst'] = hurst
#write the correlation in the table
for i in range (0, len(hurst_df)):
tbl.loc[(tbl['sub'] == sub) & (tbl['run'] == sess) & (tbl['electrodes'] == hurst_df.iloc[i, 0]), 'Hurst'] = hurst_df.iloc[i, 1]
#save table
tbl.to_excel(pathtbl, index=False)
# Plot the results
matplotlib.rcParams.update({'font.size': 11})
plt.figure(figsize=(20,5),layout='constrained')
plot_tetrode(hurst,chRegions)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_title('Hurst exponent - Memory of time series',fontsize=18)
ax.set_ylabel('H values',fontsize=15)
plt.xticks(rotation = 45)
if save==1:
plt.savefig(path+'Hurst_Exponent_AllChannels.jpg')