Application of Machine Learning Models to Identify Differences in Neural Electrophysiological Properties Across Estrous Cycle Phases.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches

[1] § Materials and Methods › Data Preprocessing and Feature Extraction: Passive Properties ↔ Estrous_Classification.ipynb, lines 657–658 · score 0.94 · voltage deflection onset, inward rectification ratio, steady state voltage, rebound voltage, sag amplitude, voltage change
[2] § Materials and Methods › Model Training and Classification ↔ Estrous_Classification.ipynb, lines 702–751 · score 0.91 · gradient boosting classifier, decision tree classifier, logistic regression, random forest classifier, neighbors classifier, score
[3] § Materials and Methods › Data Preprocessing and Feature Extraction: AP Properties ↔ Estrous_Classification.ipynb, lines 134–169 · score 0.85 · find_peaks, spike amplitude, spike peak, distance, prominence, Scipy
[4] § Materials and Methods › Data Preprocessing and Feature Extraction: mEPSC Properties ↔ Estrous_Classification.ipynb, lines 590–591 · score 0.81 · 10–90, inter event, half width, slope, noise, decay

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 2,588 lines · 88 KB · MIT · 4 matches

# %% [markdown]
# <a href="https://colab.research.google.com/github/Armaan-Raina/Estrous-Phase-Classification/blob/main/Estrous_Classification.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
# %% [markdown]
# # **Estrous Cycle Phase Classification**
# %% [markdown]
# ## **Preprocessing/ Feature Extraction**
# %% [markdown]
# ### **Importing Packages**
# %%
!pip install pyABF efel
# %%
import pyabf
import numpy
import pandas
import matplotlib.pyplot
import os
import glob
import scipy.signal
import scipy
import seaborn
import efel
import sklearn.manifold.TSNE
import sklearn.preprocessing.StandardScaler
import sklearn.model_selection.train_test_split
import sklearn.ensemble.RandomForestClassifier
import sklearn.ensemble.GradientBoostingClassifier
import sklearn.linear_model.LogisticRegression
import sklearn.neighbors.KNeighborsClassifier
import sklearn.neural_network.MLPClassifier
import sklearn.svm.SVC
import sklearn.tree.DecisionTreeClassifier
import sklearn.decomposition.PCA
import sklearn.metrics.accuracy_score
import sklearn.metrics.classification_report
import statsmodels.formula.api
import statsmodels.api
import matplotlib.colors
import pickle
import statsmodels.stats.multicomp.pairwise_tukeyhsd
import networkx as nx
import itertools
from PIL import Image
# %%
!pip freeze > requirements.txt
# %% [markdown]
# ### **Uploading raw data**
# %%
raw_mEPSC_diestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Diestrus/*.abf')
raw_mEPSC_estrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Estrus/*.abf')
raw_mEPSC_Eproestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Early Proestrus/*.abf')
raw_mEPSC_Lproestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Late Proestrus/*.abf')
# %%
raw_spike_diestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Diestrus Spiking/*.abf')
raw_spike_estrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Estrus Spiking/*.abf')
raw_spike_Eproestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Early Proestrus Spiking/*.abf')
raw_spike_Lproestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Late Proestrus Spiking/*.abf')
# %%
raw_passive_diestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Diestrus Passive/*.abf')
raw_passive_estrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Estrus Passive/*.abf')
raw_passive_Eproestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Early Proestrus Passive/*.abf')
raw_passive_Lproestrus = glob.glob('/content/drive/MyDrive/Estrous Classification Project/Late Proestrus Passive/*.abf')
# %% [markdown]
# ### **Defining functions for preprocessing**
# %%
def read_abf(filename:str):
abf = pyabf.ABF(filename)
time = np.array([], dtype=float)
current = np.array([], dtype=float)
sweep_len = abf.sweepX[-1]
for sweep in range(abf.sweepCount):
abf.setSweep(sweep)
current = np.concatenate((current, abf.sweepY))
time = np.concatenate((time, abf.sweepX + sweep * sweep_len))
return current #No longer returning time
# %%
def filter_signal(current, freq=20000):
raw_signal = current
adjusted_signal = raw_signal - np.median(raw_signal)
b_lowpass, a_lowpass = signal.bessel(4, 1000, 'low', analog=False, norm='phase', fs=freq)
b_notch, a_notch = signal.iirnotch(60.0, 30.0, fs=freq)
b_multiband = signal.convolve(b_lowpass, b_notch)
a_multiband = signal.convolve(a_lowpass, a_notch)
filtered_signal = signal.filtfilt(b_multiband, a_multiband, adjusted_signal)
return filtered_signal
# %% [markdown]
# ### **Preprocessing for Spiking Data**
# %%
abf = pyabf.ABF(raw_spike_diestrus[0])
for sweep in abf.sweepList:
abf.setSweep(sweep)
plt.plot(abf.sweepX, abf.sweepY)
# %%
#abf -> list of np.ndarrays
clean_diestrus = []
clean_estrus = []
clean_Eproestrus = []
clean_Lproestrus = []
for file in raw_spike_diestrus:
clean_diestrus.append(read_abf(file))
for file in raw_spike_estrus:
clean_estrus.append(read_abf(file))
for file in raw_spike_Eproestrus:
clean_Eproestrus.append(read_abf(file))
for file in raw_spike_Lproestrus:
clean_Lproestrus.append(read_abf(file))
# %% [markdown]
# #### **Following code snippets modified based on implemented from [link text](https://)https://spikesandbursts.wordpress.com/2023/08/24/patch-clamp-data-analysis-in-python-bursts/**
# %%
def spike_detection(sig, fs):
thresh_min = -25
thresh_prominence = 15
thresh_min_width = 0.5 * (fs/1000)
distance_min = 1 * (fs/1000)
peaks, peaks_dict = scipy.signal.find_peaks(sig,
height=thresh_min,
threshold=thresh_min,
distance=distance_min,
prominence=thresh_prominence,
width=thresh_min_width,
wlen=None, # Window length to calculate prominence
rel_height=0.5, # Relative height at which the peak width is measured
plateau_size=None)
spikes_table = pd.DataFrame(columns = ['spike', 'spike_index', 'spike_time',
'inst_freq', 'isi_s',
'width', 'rise_half_ms', 'decay_half_ms',
'spike_peak', 'spike_amplitude'])
spikes_table.spike = np.arange(1, len(peaks) + 1)
spikes_table.spike_index = peaks
spikes_table.spike_time = peaks / fs # Divided by fs to get s
spikes_table.isi_s = np.diff(peaks, axis=0, prepend=peaks[0]) / fs
spikes_table.inst_freq = 1 / spikes_table.isi_s
spikes_table.width = peaks_dict['widths']/(fs/1000) # Width (ms) at half-height
spikes_table.rise_half_ms = (peaks - peaks_dict['left_ips'])/(fs/1000)
spikes_table.decay_half_ms = (peaks_dict['right_ips'] - peaks)/(fs/1000)
spikes_table.spike_peak = peaks_dict['peak_heights'] # height parameter is needed
spikes_table.spike_amplitude = peaks_dict['prominences'] # prominence parameter is needed
return spikes_table
# %%
abf = pyabf.ABF(raw_spike_diestrus[0])
print(abf.dataPointsPerMs*1000)
# %%
##Testing on one file
df = spike_detection(clean_diestrus[0][1], 20000)
# %%
mod_df = df.drop(["spike", "spike_index", "spike_time"], axis=1)
# %%
mod_df = mod_df.replace([-np.inf, np.inf], np.nan).dropna()
# %%
plt.matshow(mod_df.corr())
cb = plt.colorbar()
cb.ax.tick_params(labelsize=14)
plt.title('Correlation Matrix', fontsize=16);
# %%
#Estimating parameters based on diestrus trace
hist_data = mod_df['isi_s']
hist_stats = pd.DataFrame()
bin_size = 10 #ms
isi_range = np.ptp(hist_data)
bins = int((isi_range * 1000 / bin_size) + 0.5) # Round to the nearest integer
hist = np.histogram(hist_data, bins=bins)
hist_counts = hist[0]
hist_bins = hist[1]
# Cumulative moving average
cum = np.cumsum(hist_counts) # Cumulative sum
cma = cum / np.arange(1, len(cum) + 1)
# Calculate peaks and valleys of the cma
cma_peaks_indexes = scipy.signal.argrelextrema(cma, np.greater)
cma_valleys_indexes = scipy.signal.argrelextrema(cma, np.less)
# Select the peak you're interested in
peak_index = cma_peaks_indexes[0][0] # Change second number to select the peak
alpha = cma[peak_index] * 0.5 # Half-peak, adapt the value to your threshold criterion
# Calculate cma_threshold_index relative to the selected cma_peak
cma_threshold = (np.argmin(cma[peak_index:] >= alpha) + peak_index) * bin_size/1000
# Dataframe with histogram statistics
length = len(hist_stats)
hist_stats.loc[length, 'mean_isi'] = np.mean(hist_data)
hist_stats.loc[length, 'median_isi'] = np.median(hist_data)
hist_stats.loc[length, 'kurtosis'] = scipy.stats.kurtosis(hist_counts)
hist_stats.loc[length, 'skewness'] = scipy.stats.skew(hist_counts, bias=True)
hist_stats.loc[length, 'cma_threshold'] = cma_threshold
hist_stats.loc[length, 'cma_valley_time'] = cma_valleys_indexes[0][1] * bin_size/1000 # Change peak index as needed
hist_stats.loc[length, 'cma_peak_time'] = cma_peaks_indexes[0][0] * bin_size/1000 # Change peak index as needed
# Plot ISI histogram
fig, ax = plt.subplots(figsize=(8, 4))
ax.set_title("ISI histogram")
ax.hist(hist_data, bins=bins, alpha=0.6)
# Plot CMA
cma_x = np.linspace(np.min(hist_bins), np.max(hist_bins), bins)
ax.plot(cma_x, cma)
# Plot CMA threshold line
ax.axvline(cma_threshold, linestyle="dotted", color="gray")
# Plot CMA valleys
ax.plot(cma_x[cma_valleys_indexes], cma[cma_valleys_indexes], 'ko')
ax.plot(cma_x[cma_peaks_indexes], cma[cma_peaks_indexes], 'mo')
# ax.set_xscale('log') # Logarithmic scale may be easier to set the threshold
ax.set_xlabel("Time bins (s)")
ax.set_ylabel("Count")
ax.set_xlim([0,1])
# Show graph and table
plt.show()
hist_stats
# %%
def burst_detection(df, spike_times, spike_amplitudes, spike_peaks, n_spikes, max_isi, min_ibi=None):
df = df.sort_values(by=spike_times)
df['burst'] = np.nan
burst_num = 0
burst_start = None
last_spike = None
for i, row in df.iterrows(): # Loop through DataFrame rows
spike = row[spike_times] # Extract the spike position
if burst_start is None: # It checks if it is the first spike
burst_start = spike # It marks the current spike position as the start of a burst
last_spike = spike # Update the last_spike position to the current spike position
df.at[i, 'burst'] = burst_num # Assign burst number
elif spike - last_spike <= max_isi: # It checks if the current spike is within max isi
df.at[i, 'burst'] = burst_num
last_spike = spike
elif spike - last_spike > min_ibi: # It checks if the interburst interval has been reached
burst_num += 1
burst_start = spike
last_spike = spike
df.at[i, 'burst'] = burst_num
# Filter bursts with less than min_spikes
df = df[df.groupby('burst')[spike_times].transform('count') >= n_spikes]
bursts = df.groupby('burst')[spike_times].agg(['min', 'max', 'count'])
bursts.columns = ['burst_start', 'burst_end', 'spikes_in_bursts']
bursts['burst_length'] = bursts['burst_end'] - bursts['burst_start']
bursts['avg_spike_amplitude'] = df.groupby('burst')[spike_amplitudes].mean()
bursts['avg_spike_peaks'] = df.groupby('burst')[spike_peaks].mean()
bursts['spikes_frequency'] = bursts['spikes_in_bursts'] / bursts['burst_length']
bursts = bursts.reset_index()
bursts['burst_number'] = bursts.index + 1
return bursts[['burst_number', 'burst_start', 'burst_end',
'burst_length', 'spikes_in_bursts', 'avg_spike_amplitude',
'avg_spike_peaks', 'spikes_frequency']]
# %%
diestrus_spiking = pd.DataFrame()
estrus_spiking = pd.DataFrame()
Eproestrus_spiking = pd.DataFrame()
Lproestrus_spiking = pd.DataFrame()
for trace in clean_diestrus:
spikes_table = spike_detection(trace[1], 20000)
diestrus_spiking = pd.concat([diestrus_spiking, burst_detection(spikes_table, 'spike_time', 'spike_amplitude', 'spike_peak', n_spikes = 2, max_isi=0.5, min_ibi=0.1)])
for trace in clean_estrus:
spikes_table = spike_detection(trace[1], 20000)
estrus_spiking = pd.concat([estrus_spiking, burst_detection(spikes_table, 'spike_time', 'spike_amplitude', 'spike_peak', n_spikes = 2, max_isi=0.5, min_ibi=0.1)])
for trace in clean_Eproestrus:
spikes_table = spike_detection(trace[1], 20000)
Eproestrus_spiking = pd.concat([Eproestrus_spiking, burst_detection(spikes_table, 'spike_time', 'spike_amplitude', 'spike_peak', n_spikes = 2, max_isi=0.5, min_ibi=0.1)])
for trace in clean_Lproestrus:
spikes_table = spike_detection(trace[1], 20000)
Lproestrus_spiking = pd.concat([Lproestrus_spiking, burst_detection(spikes_table, 'spike_time', 'spike_amplitude', 'spike_peak', n_spikes = 2, max_isi=0.5, min_ibi=0.1)])
# %%
diestrus_spiking.to_csv("/content/drive/MyDrive/Estrous Classification Project/diestrus_spiking_feats")
estrus_spiking.to_csv("/content/drive/MyDrive/Estrous Classification Project/estrus_spiking_feats")
Eproestrus_spiking.to_csv("/content/drive/MyDrive/Estrous Classification Project/Earlyproestrus_spiking_feats")
Lproestrus_spiking.to_csv("/content/drive/MyDrive/Estrous Classification Project/Lateproestrus_spiking_feats")
# %% [markdown]
# ### **Preprocessing for Passive Data**
# %%
#Exploration of abf file to determine constants
abf = pyabf.ABF(raw_passive_Lproestrus[0])
sweeps = abf.sweepList
for sweep in sweeps:
abf.setSweep(sweep)
plt.plot(abf.sweepX, abf.sweepC)
abf.setSweep(4)
##Finding the start and end indices of the stimulation
stim_start = 0
stim_end = 0
idx = 0
startFlag = False
val = abf.sweepC
while val[idx] == 0:
idx+=1
stim_start += 1
print(stim_start)
for i in range(stim_start, len(val)+1):
if val[i] == 0:
stim_end = i
break
print(stim_end)
print(abf.sweepX[9562])
print(abf.sweepX[stim_end])
plt.plot(abf.sweepX, abf.sweepC)
plt.axvline(0.4781, color='r')
plt.axvline(1.0781, color='g')
print(abf.sweepEpochs)
# %%
##CONSTANTS - determined in last block
stim_start = [4781 / abf.dataRate * 1000] # ms
stim_end = [10781 / abf.dataRate * 1000] # ms
ir_current_start = [i for i in range(-140, 10, 10)]
ir_current_stop = 0
curr_start_idx = 4781
curr_stop_idx = 10781
# %%
ir_current_start, abf.dataPointsPerMs * 1000 #20 kHz sampling rate
# %%
def init_features(filename):
abf = pyabf.ABF(filename)
sweeps = abf.sweepList
table = pd.DataFrame(columns=[
'rmp_mV',
'steady_state_voltage_mV',
'vmin_mV',
'vdeflection_begin_mV',
'vrebound_mV',
'voltage_delta_mV',
'current_pA',
'input_resistance_Gohm',
'time_constant_ms',
'capacitance_pF',
'sag_amplitude_mV',
'sag_ratio1',
'sag_ratio2',
'sweep',
'inward_rectification_ratio'
])
inwardRec = []
inwardRecRatio = np.nan
for sweep in sweeps:
abf.setSweep(sweep)
stim_start_index = 4781
stim_end_index = 10781
stim_start = [abf.sweepX[stim_start_index] * 1000]
stim_end = [abf.sweepX[stim_end_index] * 1000]
current_pA = np.mean(abf.sweepC[stim_start_index:stim_end_index])
trace = {
'T': abf.sweepX * 1000,
'V': abf.sweepY,
'stim_start': stim_start,
'stim_end': stim_end
}
feature_values = efel.get_feature_values(
[trace],
[
'voltage_base', 'steady_state_voltage_stimend',
'minimum_voltage', 'voltage_deflection_begin',
'voltage_deflection', 'voltage_deflection_vb_ssse',
'decay_time_constant_after_stim', 'sag_amplitude',
'sag_ratio1', 'sag_ratio2', 'voltage_after_stim'
]
)[0]
def safe_get(key, idx=0): #preventing nulls so it can be handled later
return feature_values[key][idx] if (key in feature_values and feature_values[key] is not None) else np.nan
base_v = safe_get('voltage_base')
steady_v = safe_get('steady_state_voltage_stimend')
vmin = safe_get('minimum_voltage')
vdeflection_begin = safe_get('voltage_deflection_begin')
vrebound = safe_get('voltage_after_stim')
sag_amp = safe_get('sag_amplitude')
sag_ratio1 = safe_get('sag_ratio1')
sag_ratio2 = safe_get('sag_ratio2')
tau = safe_get('decay_time_constant_after_stim')
delta_v = steady_v - base_v if not np.isnan(steady_v) and not np.isnan(base_v) else np.nan
Rin = delta_v / current_pA if current_pA != 0 else np.nan
capacitance = tau / Rin if Rin and not np.isnan(tau) and Rin != 0 else np.nan
# record specific sweeps for inward rectification ratio
if sweep == 1:
inwardRec.append(Rin)
if sweep == 14:
inwardRec.append(Rin)
table.loc[len(table)] = [
base_v, steady_v, vmin, vdeflection_begin, vrebound, delta_v,
current_pA, Rin, tau, capacitance, sag_amp, sag_ratio1, sag_ratio2, sweep, np.nan
]
if len(inwardRec) == 2 and not any(np.isnan(inwardRec)):
inwardRecRatio = inwardRec[1] / inwardRec[0]
table.iloc[-len(sweeps):, table.columns.get_loc('inward_rectification_ratio')] = inwardRecRatio
return table
# %%
diestrus_passive_df = pd.DataFrame()
estrus_passive_df = pd.DataFrame()
Eproestrus_passive_df = pd.DataFrame()
Lproestrus_passive_df = pd.DataFrame()
for filename in raw_passive_diestrus:
df = init_features(filename)
diestrus_passive_df = pd.concat([df, diestrus_passive_df])
diestrus_passive_df['label'] = 0
for filename in raw_passive_estrus:
df = init_features(filename)
estrus_passive_df = pd.concat([df, estrus_passive_df])
estrus_passive_df['label'] = 1
for filename in raw_passive_Eproestrus:
df = init_features(filename)
Eproestrus_passive_df = pd.concat([df, Eproestrus_passive_df])
Eproestrus_passive_df['label'] = 2
for filename in raw_passive_Lproestrus:
df = init_features(filename)
Lproestrus_passive_df = pd.concat([df, Lproestrus_passive_df])
Lproestrus_passive_df['label'] = 3
# %%
diestrus_passive_df.head(n = 45)
# %%
diestrus_passive_df = diestrus_passive_df.dropna()
estrus_passive_df = estrus_passive_df.dropna()
Eproestrus_passive_df = Eproestrus_passive_df.dropna()
Lproestrus_passive_df = Lproestrus_passive_df.dropna()
# %%
diestrus_passive_df.to_csv("/content/drive/MyDrive/Estrous Classification Project/diestrus_passive_feats")
estrus_passive_df.to_csv("/content/drive/MyDrive/Estrous Classification Project/estrus_passive_feats")
Eproestrus_passive_df.to_csv("/content/drive/MyDrive/Estrous Classification Project/Eproestrus_passive_feats")
Lproestrus_passive_df.to_csv("/content/drive/MyDrive/Estrous Classification Project/Lproestrus_passive_feats")
# %% [markdown]
# ## **Feature Loading**
#
# %% [markdown]
# ### **mEPSC (MiniAnalysis) Features**
# %%
mEPSC_diestrus = pd.read_excel('/content/drive/MyDrive/Spreadsheets/diestrus.xlsx')
mEPSC_estrus = pd.read_excel('/content/drive/MyDrive/Spreadsheets/estrus.xlsx')
mEPSC_Eproestrus = pd.read_excel('/content/drive/MyDrive/Spreadsheets/earlyproestrus.xlsx')
mEPSC_Lproestrus = pd.read_excel('/content/drive/MyDrive/Spreadsheets/lateproestrus.xlsx')
# %%
#Filtering out extraneous events likely due to experimental setup or other circumstances
mEPSC_diestrus = mEPSC_diestrus[(mEPSC_diestrus['Amplitude'] < 50) & (mEPSC_diestrus['Baseline'] > 0)]
mEPSC_estrus = mEPSC_estrus[(mEPSC_estrus['Amplitude'] < 50) & (mEPSC_estrus['Baseline'] > 0)]
mEPSC_Eproestrus = mEPSC_Eproestrus[(mEPSC_Eproestrus['Amplitude'] < 50) & (mEPSC_Eproestrus['Baseline'] > 0)]
mEPSC_Lproestrus = mEPSC_Lproestrus[(mEPSC_Lproestrus['Amplitude'] < 50) & (mEPSC_Lproestrus['Baseline'] > 0)]
# %%
def compute_iei(df: pd.DataFrame) -> pd.DataFrame:
df["Inter-Event Interval (ms)"] = df["Time (ms)"].diff()
return df
# %%
#Simple IEI computation
mEPSC_diestrus = compute_iei(mEPSC_diestrus)
mEPSC_estrus = compute_iei(mEPSC_estrus)
mEPSC_Eproestrus = compute_iei(mEPSC_Eproestrus)
mEPSC_Lproestrus = compute_iei(mEPSC_Lproestrus)
# %%
mEPSC_diestrus.iloc[0, 18] = mEPSC_diestrus.iloc[0, 1]
mEPSC_estrus.iloc[0, 18] = mEPSC_estrus.iloc[0, 1]
mEPSC_Eproestrus.iloc[0, 18] = mEPSC_Eproestrus.iloc[0, 1]
mEPSC_Lproestrus.iloc[0, 18] = mEPSC_Lproestrus.iloc[0, 1]
# %%
mEPSC_min_count = len(min([mEPSC_diestrus, mEPSC_estrus, mEPSC_Eproestrus, mEPSC_Lproestrus], key=len))
# %%
mEPSC_diestrus['label'] = 0
mEPSC_estrus['label'] = 1
mEPSC_Eproestrus['label'] = 2
mEPSC_Lproestrus['label'] = 3
# %%
mEPSC_diestrus = mEPSC_diestrus.dropna()
mEPSC_estrus = mEPSC_estrus.dropna()
mEPSC_Eproestrus = mEPSC_Eproestrus.dropna()
mEPSC_Lproestrus = mEPSC_Lproestrus.dropna()
# %%
mEPSC = pd.concat([mEPSC_diestrus, mEPSC_estrus, mEPSC_Eproestrus, mEPSC_Lproestrus])
# %%
#Reducing count of each class to match for testing to determine if class imbalance significantly impacts model accuracy
mEPSC_even = pd.concat([mEPSC_diestrus.head(mEPSC_min_count), mEPSC_estrus.head(mEPSC_min_count), mEPSC_Eproestrus.head(mEPSC_min_count), mEPSC_Lproestrus.head(mEPSC_min_count)])
# %%
mEPSC_even.shape
# %%
mEPSC_even = mEPSC_even.drop(["Obs Num", "Time (ms)", "Group", "Channel", "Peak Dir", "Burst#", "BurstE#", "Rel Time"], axis=1)
# %%
mEPSC = mEPSC.drop(["Obs Num", "Time (ms)", "Group", "Channel", "Peak Dir", "Burst#", "BurstE#", "Rel Time"], axis=1) #Dropping insignificant features
# %%
mEPSC_names = ["Amplitude", "Rise (ms)", "Decay (ms)", "Area", "Baseline", "Noise", "10-90Rise", "HalfWidth", "Rise50", "10-90Slope", "Inter-Event Interval (ms)"]
# %% [markdown]
# ### **AP Features**
# %%
spiking_diestrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/diestrus_spiking_feats")
spiking_estrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/estrus_spiking_feats")
spiking_Eproestrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/Earlyproestrus_spiking_feats")
spiking_Lproestrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/Lateproestrus_spiking_feats")
# %%
spiking_min_count = len(min([spiking_diestrus, spiking_estrus, spiking_Eproestrus, spiking_Lproestrus], key=len))
# %%
spiking_diestrus.shape, spiking_estrus.shape, spiking_Eproestrus.shape, spiking_Lproestrus.shape
# %%
spiking_diestrus['label'] = 0
spiking_estrus['label'] = 1
spiking_Eproestrus['label'] = 2
spiking_Lproestrus['label'] = 3
# %%
spiking = pd.concat([spiking_diestrus, spiking_estrus, spiking_Eproestrus, spiking_Lproestrus])
# %%
spiking = spiking.drop(["Unnamed: 0", "burst_number", "burst_start", "burst_end"], axis=1) #Dropping insignificant features
# %%
spiking_even = pd.concat([spiking_diestrus.head(spiking_min_count), spiking_estrus.head(spiking_min_count), spiking_Eproestrus.head(spiking_min_count), spiking_Lproestrus.head(spiking_min_count)])
# %%
spiking_even = spiking_even.drop(["Unnamed: 0", "burst_number", "burst_start", "burst_end"], axis=1) #Dropping insignificant features
# %%
spiking_names = ["burst_length", "spikes_in_bursts", "avg_spike_amplitude", "avg_spike_peaks", "spikes_frequency"]
# %% [markdown]
# ### **Passive Features**
# %%
passive_diestrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/diestrus_passive_feats")
passive_estrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/estrus_passive_feats")
passive_Eproestrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/Eproestrus_passive_feats")
passive_Lproestrus = pd.read_csv("/content/drive/MyDrive/Estrous Classification Project/Lproestrus_passive_feats")
#These already have labels saved into them
# %%
passive_min_count = len(min([passive_diestrus, passive_estrus, passive_Eproestrus, passive_Lproestrus], key=len))
# %%
passive = pd.concat([passive_diestrus, passive_estrus, passive_Eproestrus, passive_Lproestrus])
# %%
passive = passive.drop(["Unnamed: 0", "sag_ratio1", "sweep"], axis=1)
# %%
passive_even = pd.concat([passive_diestrus.head(passive_min_count), passive_estrus.head(passive_min_count), passive_Eproestrus.head(passive_min_count), passive_Lproestrus.head(passive_min_count)])
# %%
passive_even = passive_even.drop(["Unnamed: 0", "sag_ratio1", "sweep"], axis=1)
# %%
passive_names = ["rmp_mV", "steady_state_voltage_mV", "vmin_mV", "vdeflection_begin_mV", "vrebound_mV", "voltage_delta_mV", "current_pA", "input_resistance_Gohm", "time_constant_ms", "capacitance_pF", "sag_amplitude_mV", "sag_ratio1", "sag_ratio2", "sweep", "inward_rectification_ratio"]
# %%
passive_names_pretty = ["RMP (mV)", "Steady State Voltage (mV)", "Min. Voltage (mV)", "Voltage Deflection Onset (mV)", "Rebound Voltage (mV)", "Voltage Change (mV)", "Current (pA)", "Input Resistance (GOhms)", "time_constant_ms", "capacitance_pF", "sag_amplitude_mV","sag_ratio2", "inward_rectification_ratio"]
# %% [markdown]
# ## **Model Training**
# %% [markdown]
# #### **Feature Preparation**
# %%
X_mEPSC = mEPSC.drop('label', axis=1)
y_mEPSC = mEPSC['label']
# %%
X_spiking = spiking.drop('label', axis=1)
y_spiking = spiking['label']
# %%
X_passive = passive.drop('label', axis=1)
y_passive = passive['label']
# %%
X_even_mEPSC = mEPSC_even.drop('label', axis=1)
y_even_mEPSC = mEPSC_even['label']
X_even_spiking = spiking_even.drop('label', axis=1)
y_even_spiking = spiking_even['label']
X_even_passive = passive_even.drop('label', axis=1)
y_even_passive = passive_even['label']
# %%
sc = StandardScaler()
X_mEPSC = sc.fit_transform(X_mEPSC)
X_spiking = sc.fit_transform(X_spiking)
X_passive = sc.fit_transform(X_passive)
X_even_mEPSC = sc.fit_transform(X_even_mEPSC)
X_even_spiking = sc.fit_transform(X_even_spiking)
X_even_passive = sc.fit_transform(X_even_passive)
# %% [markdown]
# #### **mEPSC First**
# %%
mEPSC_classification_reports = {}
mEPSC_feature_importances = {}
# %%
SEEDS = [42, 7, 13, 21, 99]
mEPSC_seed_accuracies = {model_name: [] for model_name in [
'RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC'
]}
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
# Store ALL predictions and true labels for ALL seeds' runs
all_preds_mEPSC = {name: [] for name in model_names}
all_y_tests_mEPSC = {name: [] for name in model_names}
for seed in SEEDS:
X_train_mEPSC, X_test_mEPSC, y_train_mEPSC, y_test_mEPSC = train_test_split(
X_mEPSC, y_mEPSC, test_size=0.25, shuffle=True, random_state=seed
)
mEPSC_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(), # No random_state
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, mEPSC_models):
model.fit(X_train_mEPSC, y_train_mEPSC)
pred = model.predict(X_test_mEPSC)
acc = accuracy_score(y_test_mEPSC, pred)
mEPSC_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_mEPSC[name].extend(pred)
all_y_tests_mEPSC[name].extend(y_test_mEPSC)
# Feature importances for each seed of the RFC
if name == 'RFC':
mEPSC_feature_importances[seed] = model.feature_importances_
# Average across seeds
mEPSC_accuracies = [np.mean(mEPSC_seed_accuracies[name]) for name in model_names]
mEPSC_std = [np.std(mEPSC_seed_accuracies[name]) for name in model_names]
for name in model_names:
if name in all_preds_mEPSC and name in all_y_tests_mEPSC:
rep = classification_report(all_y_tests_mEPSC[name], all_preds_mEPSC[name], zero_division=0, output_dict=True)
mEPSC_classification_reports[name] = rep
# %% [markdown]
# #### **Spiking Next**
# %%
spiking_classification_reports = {}
spiking_feature_importances = {}
# %%
SEEDS = [42, 7, 13, 21, 99] # Same seeds
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
spiking_seed_accuracies = {name: [] for name in model_names}
all_preds_spiking = {name: [] for name in model_names}
all_y_tests_spiking = {name: [] for name in model_names}
for seed in SEEDS:
X_train_spiking, X_test_spiking, y_train_spiking, y_test_spiking = train_test_split(
X_spiking, y_spiking, test_size=0.25, shuffle=True, random_state=seed
)
spiking_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(), # No random_state
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, spiking_models):
model.fit(X_train_spiking, y_train_spiking)
pred = model.predict(X_test_spiking)
acc = accuracy_score(y_test_spiking, pred)
spiking_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_spiking[name].extend(pred)
all_y_tests_spiking[name].extend(y_test_spiking)
if name == 'RFC':
spiking_feature_importances[seed] = model.feature_importances_
# Average across seeds
spiking_accuracies = [np.mean(spiking_seed_accuracies[name]) for name in model_names]
spiking_std = [np.std(spiking_seed_accuracies[name]) for name in model_names]
for name in model_names:
if name in all_preds_spiking and name in all_y_tests_spiking:
rep = classification_report(all_y_tests_spiking[name], all_preds_spiking[name], zero_division=0, output_dict=True)
spiking_classification_reports[name] = rep
# %% [markdown]
# #### **Passive Last**
# %%
passive_classification_reports = {}
passive_feature_importances = {}
# %%
SEEDS = [42, 7, 13, 21, 99]
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
passive_seed_accuracies = {name: [] for name in model_names}
# Store ALL predictions and true labels for ALL seeds' runs
all_preds_passive = {name: [] for name in model_names}
all_y_tests_passive = {name: [] for name in model_names}
for seed in SEEDS:
X_train_passive, X_test_passive, y_train_passive, y_test_passive = train_test_split(
X_passive, y_passive, test_size=0.25, shuffle=True, random_state=seed
)
passive_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(),
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, passive_models):
model.fit(X_train_passive, y_train_passive)
pred = model.predict(X_test_passive)
acc = accuracy_score(y_test_passive, pred)
passive_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_passive[name].extend(pred)
all_y_tests_passive[name].extend(y_test_passive)
if name == "RFC":
passive_feature_importances[seed] = model.feature_importances_
# Average across seeds
passive_accuracies = [np.mean(passive_seed_accuracies[name]) for name in model_names]
passive_std = [np.std(passive_seed_accuracies[name]) for name in model_names]
for name in model_names:
if name in all_preds_passive and name in all_y_tests_passive:
rep = classification_report(all_y_tests_passive[name], all_preds_passive[name], zero_division=0, output_dict=True)
passive_classification_reports[name] = rep
# %% [markdown]
# #### **Now again, addressing data inequality**
# %%
###mEPSC###
SEEDS = [42, 7, 13, 21, 99]
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
mEPSC_even_seed_accuracies = {name: [] for name in model_names}
all_preds_mEPSC_even = {name: [] for name in model_names}
all_y_tests_mEPSC_even = {name: [] for name in model_names}
for seed in SEEDS:
X_train_mEPSC_even, X_test_mEPSC_even, y_train_mEPSC_even, y_test_mEPSC_even = train_test_split(
X_even_mEPSC, y_even_mEPSC, test_size=0.25, shuffle=True, random_state=seed
)
mEPSC_even_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(),
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, mEPSC_even_models):
model.fit(X_train_mEPSC_even, y_train_mEPSC_even)
pred = model.predict(X_test_mEPSC_even)
acc = accuracy_score(y_test_mEPSC_even, pred)
mEPSC_even_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_mEPSC_even[name].extend(pred)
all_y_tests_mEPSC_even[name].extend(y_test_mEPSC_even)
# Average across seeds
mEPSC_even_accuracies = [np.mean(mEPSC_even_seed_accuracies[name]) for name in model_names]
mEPSC_even_std = [np.std(mEPSC_even_seed_accuracies[name]) for name in model_names]
# %%
###Spiking###
SEEDS = [42, 7, 13, 21, 99]
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
spiking_even_seed_accuracies = {name: [] for name in model_names}
# Store ALL predictions and true labels for ALL seeds' runs
all_preds_spiking_even = {name: [] for name in model_names}
all_y_tests_spiking_even = {name: [] for name in model_names}
for seed in SEEDS:
X_train_spiking_even, X_test_spiking_even, y_train_spiking_even, y_test_spiking_even = train_test_split(
X_even_spiking, y_even_spiking, test_size=0.25, shuffle=True, random_state=seed
)
spiking_even_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(),
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, spiking_even_models):
model.fit(X_train_spiking_even, y_train_spiking_even)
pred = model.predict(X_test_spiking_even)
acc = accuracy_score(y_test_spiking_even, pred)
spiking_even_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_spiking_even[name].extend(pred)
all_y_tests_spiking_even[name].extend(y_test_spiking_even)
# Average across seeds
spiking_even_accuracies = [np.mean(spiking_even_seed_accuracies[name]) for name in model_names]
spiking_even_std = [np.std(spiking_even_seed_accuracies[name]) for name in model_names]
# %%
###Passive###
SEEDS = [42, 7, 13, 21, 99]
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
passive_even_seed_accuracies = {name: [] for name in model_names}
all_preds_passive_even = {name: [] for name in model_names}
all_y_tests_passive_even = {name: [] for name in model_names}
for seed in SEEDS:
X_train_passive_even, X_test_passive_even, y_train_passive_even, y_test_passive_even = train_test_split(
X_even_passive, y_even_passive, test_size=0.25, shuffle=True, random_state=seed
)
passive_even_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(),
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, passive_even_models):
model.fit(X_train_passive_even, y_train_passive_even)
pred = model.predict(X_test_passive_even)
acc = accuracy_score(y_test_passive_even, pred)
passive_even_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_passive_even[name].extend(pred)
all_y_tests_passive_even[name].extend(y_test_passive_even)
# Average across seeds
passive_even_accuracies = [np.mean(passive_even_seed_accuracies[name]) for name in model_names]
passive_even_std = [np.std(passive_even_seed_accuracies[name]) for name in model_names]
# %% [markdown]
# #### **Once more with randomized labels**
# %%
rng = np.random.default_rng(seed=42)
# %%
SEEDS = [42, 7, 13, 21, 99]
mEPSC_rnd_seed_accuracies = {model_name: [] for model_name in [
'RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC'
]}
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
all_preds_mEPSC_rnd = {name: [] for name in model_names}
all_y_tests_mEPSC_rnd = {name: [] for name in model_names}
for seed in SEEDS:
X_train_mEPSC, X_test_mEPSC, y_train_mEPSC, y_test_mEPSC = train_test_split(
X_mEPSC, y_mEPSC, test_size=0.25, shuffle=True, random_state=seed
)
y_train_mEPSC_rnd = rng.permutation(y_train_mEPSC)
mEPSC_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(), # No random_state
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, mEPSC_models):
model.fit(X_train_mEPSC, y_train_mEPSC_rnd)
pred = model.predict(X_test_mEPSC)
acc = accuracy_score(y_test_mEPSC, pred)
mEPSC_rnd_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_mEPSC_rnd[name].extend(pred)
all_y_tests_mEPSC_rnd[name].extend(y_test_mEPSC)
# Average across seeds (accuracy)
mEPSC_rnd_accuracies = [np.mean(mEPSC_rnd_seed_accuracies[name]) for name in model_names]
mEPSC_rnd_std = [np.std(mEPSC_rnd_seed_accuracies[name]) for name in model_names]
# %%
SEEDS = [42, 7, 13, 21, 99] # Same seeds as mEPSC for consistency
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
spiking_rnd_seed_accuracies = {name: [] for name in model_names}
all_preds_spiking_rnd = {name: [] for name in model_names}
all_y_tests_spiking_rnd = {name: [] for name in model_names}
for seed in SEEDS:
X_train_spiking, X_test_spiking, y_train_spiking, y_test_spiking = train_test_split(
X_spiking, y_spiking, test_size=0.25, shuffle=True, random_state=seed
)
y_train_spiking_rnd = rng.permutation(y_train_spiking)
spiking_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(), # No random_state
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, spiking_models):
model.fit(X_train_spiking, y_train_spiking_rnd)
pred = model.predict(X_test_spiking)
acc = accuracy_score(y_test_spiking, pred)
spiking_rnd_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_spiking_rnd[name].extend(pred)
all_y_tests_spiking_rnd[name].extend(y_test_spiking)
# Average across seeds
spiking_rnd_accuracies = [np.mean(spiking_rnd_seed_accuracies[name]) for name in model_names]
spiking_rnd_std = [np.std(spiking_rnd_seed_accuracies[name]) for name in model_names]
# %%
SEEDS = [42, 7, 13, 21, 99]
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
passive_rnd_seed_accuracies = {name: [] for name in model_names}
# Store ALL predictions and true labels for ALL seeds' runs
all_preds_passive_rnd = {name: [] for name in model_names}
all_y_tests_passive_rnd = {name: [] for name in model_names}
for seed in SEEDS:
X_train_passive, X_test_passive, y_train_passive, y_test_passive = train_test_split(
X_passive, y_passive, test_size=0.25, shuffle=True, random_state=seed
)
y_train_passive_rnd = rng.permutation(y_train_passive)
passive_models = [
RandomForestClassifier(random_state=seed),
GradientBoostingClassifier(random_state=seed),
LogisticRegression(max_iter=1000, random_state=seed),
KNeighborsClassifier(),
MLPClassifier(max_iter=1000, random_state=seed),
SVC(random_state=seed),
DecisionTreeClassifier(random_state=seed),
]
for name, model in zip(model_names, passive_models):
model.fit(X_train_passive, y_train_passive_rnd)
pred = model.predict(X_test_passive)
acc = accuracy_score(y_test_passive, pred)
passive_rnd_seed_accuracies[name].append(acc)
# Store for all seeds
all_preds_passive_rnd[name].extend(pred)
all_y_tests_passive_rnd[name].extend(y_test_passive)
# Average across seeds
passive_rnd_accuracies = [np.mean(passive_rnd_seed_accuracies[name]) for name in model_names]
passive_rnd_std = [np.std(passive_rnd_seed_accuracies[name]) for name in model_names]
# %% [markdown]
# #### **Saving Results**
# %%
mEPSC_results = {
"Accuracies": mEPSC_seed_accuracies,
"Classification Reports": mEPSC_classification_reports,
"Feature Importances": mEPSC_feature_importances
}
spiking_results = {
"Accuracies": spiking_seed_accuracies,
"Classification Reports": spiking_classification_reports,
"Feature Importances": spiking_feature_importances
}
passive_results = {
"Accuracies": passive_seed_accuracies,
"Classification Reports": passive_classification_reports,
"Feature Importances": passive_feature_importances
}
mEPSC_even_results = {
"Accuracies": mEPSC_even_seed_accuracies
}
spiking_even_results = {
"Accuracies": spiking_even_seed_accuracies
}
passive_even_results = {
"Accuracies": passive_even_seed_accuracies
}
mEPSC_rnd_results = {
"Accuracies": mEPSC_rnd_seed_accuracies
}
spiking_rnd_results = {
"Accuracies": spiking_rnd_seed_accuracies
}
passive_rnd_results = {
"Accuracies": passive_rnd_seed_accuracies
}
# Define the base path for saving to Google Drive
base_path = "/content/drive/MyDrive/Estrous Classification Project/"
with open(base_path + "mEPSC_results.pkl", "wb") as f:
pickle.dump(mEPSC_results, f)
with open(base_path + "spiking_results.pkl", "wb") as f:
pickle.dump(spiking_results, f)
with open(base_path + "passive_results.pkl", "wb") as f:
pickle.dump(passive_results, f)
with open(base_path + "mEPSC_even_results.pkl", "wb") as f:
pickle.dump(mEPSC_even_results, f)
with open(base_path + "spiking_even_results.pkl", "wb") as f:
pickle.dump(spiking_even_results, f)
with open(base_path + "passive_even_results.pkl", "wb") as f:
pickle.dump(passive_even_results, f)
with open(base_path + "mEPSC_rnd_results.pkl", "wb") as f:
pickle.dump(mEPSC_rnd_results, f)
with open(base_path + "spiking_rnd_results.pkl", "wb") as f:
pickle.dump(spiking_rnd_results, f)
with open(base_path + "passive_rnd_results.pkl", "wb") as f:
pickle.dump(passive_rnd_results, f)
# %% [markdown]
# ## **Results/Analysis**
# %% [markdown]
# ### **Loading Model Training Results**
# %%
# Define the base path for loading from Google Drive
base_path = "/content/drive/MyDrive/Estrous Classification Project/"
with open(base_path + "mEPSC_results.pkl", "rb") as f:
mEPSC_results = pickle.load(f)
with open(base_path + "spiking_results.pkl", "rb") as f:
spiking_results = pickle.load(f)
with open(base_path + "passive_results.pkl", "rb") as f:
passive_results = pickle.load(f)
with open(base_path + "mEPSC_even_results.pkl", "rb") as f:
mEPSC_even_results = pickle.load(f)
with open(base_path + "spiking_even_results.pkl", "rb") as f:
spiking_even_results = pickle.load(f)
with open(base_path + "passive_even_results.pkl", "rb") as f:
passive_even_results = pickle.load(f)
with open(base_path + "mEPSC_rnd_results.pkl", "rb") as f:
mEPSC_rnd_results = pickle.load(f)
with open(base_path + "spiking_rnd_results.pkl", "rb") as f:
spiking_rnd_results = pickle.load(f)
with open(base_path + "passive_rnd_results.pkl", "rb") as f:
passive_rnd_results = pickle.load(f)
# %%
mEPSC_seed_accuracies = mEPSC_results["Accuracies"]
mEPSC_classification_reports = mEPSC_results["Classification Reports"]
mEPSC_feature_importances = mEPSC_results["Feature Importances"]
spiking_seed_accuracies = spiking_results["Accuracies"]
spiking_classification_reports = spiking_results["Classification Reports"]
spiking_feature_importances = spiking_results["Feature Importances"]
passive_seed_accuracies = passive_results["Accuracies"]
passive_classification_reports = passive_results["Classification Reports"]
passive_feature_importances = passive_results["Feature Importances"]
mEPSC_even_seed_accuracies = mEPSC_even_results["Accuracies"]
spiking_even_seed_accuracies = spiking_even_results["Accuracies"]
passive_even_seed_accuracies = passive_even_results["Accuracies"]
mEPSC_rnd_seed_accuracies = mEPSC_rnd_results["Accuracies"]
spiking_rnd_seed_accuracies = spiking_rnd_results["Accuracies"]
passive_rnd_seed_accuracies = passive_rnd_results["Accuracies"]
# %% [markdown]
# ### **Feature Importances**
#
# %%
seeds = mEPSC_feature_importances.keys()
mEPSC_imp = np.vstack([mEPSC_feature_importances[seed] for seed in seeds])
spiking_imp = np.vstack([spiking_feature_importances[seed] for seed in seeds])
passive_imp = np.vstack([passive_feature_importances[seed] for seed in seeds])
# %%
mEPSC_summary = pd.DataFrame({
"mean": mEPSC_imp.mean(axis=0),
"std": mEPSC_imp.std(axis=0),
"top5_freq": (
pd.DataFrame(mEPSC_imp, columns=mEPSC_names)
.rank(axis=1, ascending=False) <= 5
).mean(axis=0).values,
}, index=mEPSC_names).sort_values("mean", ascending=False)
passive_summary = pd.DataFrame({
"mean": passive_imp.mean(axis=0),
"std": passive_imp.std(axis=0),
"top5_freq": (
pd.DataFrame(passive_imp, columns=passive_names_pretty)
.rank(axis=1, ascending=False) <= 5
).mean(axis=0).values,
}, index=passive_names_pretty).sort_values("mean", ascending=False)
spiking_summary = pd.DataFrame({
"mean": spiking_imp.mean(axis=0),
"std": spiking_imp.std(axis=0),
"top5_freq": (
pd.DataFrame(spiking_imp, columns=spiking_names)
.rank(axis=1, ascending=False) <= 5
).mean(axis=0).values,
}, index=spiking_names).sort_values("mean", ascending=False)
# %%
mEPSC_importances = mEPSC_summary.head(5)
spiking_importances = spiking_summary.head(5)
passive_importances = passive_summary.head(5)
mEPSC_top5 = mEPSC_importances.index.tolist()
spiking_top5 = spiking_importances.index.tolist()
passive_top5 = passive_importances.index.tolist()
# %%
passive_top5 = ['Input Resistance',
'Membrane Time Constant (ms)',
'Rebound Voltage (mV)',
'Voltage Deflection Onset (mV)',
'RMP (mV)']
# %%
spiking_top5 = ['Burst Length (ms)',
'Spike Frequency',
'Average Spike Amplitude',
'Average Spike Peaks',
'Spikes per Burst']
# %%
mEPSC_top5 = ['Baseline',
'Amplitude',
'HalfWidth',
'10-90 Slope',
'Inter-Event Interval (ms)']
# %%
print("\n" + "="*60)
print("ALL FEATURES BY IMPORTANCE")
print("="*60)
print("\nmEPSC Dataset:")
print(f" {'Rank':<6} {'Feature':<35} {'Importance'}")
print(f" {'-'*6} {'-'*35} {'-'*10}")
for rank, (name, row) in enumerate(mEPSC_summary.iterrows(), start=1):
marker = " " if rank > 5 else "* "
print(f" {marker}{rank:<6} {name:<35} {row['mean']:.4f}")
print("\nSpiking Dataset:")
print(f" {'Rank':<6} {'Feature':<35} {'Importance'}")
print(f" {'-'*6} {'-'*35} {'-'*10}")
for rank, (name, row) in enumerate(spiking_summary.iterrows(), start=1):
marker = " " if rank > 5 else "* "
print(f" {marker}{rank:<6} {name:<35} {row['mean']:.4f}")
print("\nPassive Dataset:")
print(f" {'Rank':<6} {'Feature':<35} {'Importance'}")
print(f" {'-'*6} {'-'*35} {'-'*10}")
for rank, (name, row) in enumerate(passive_summary.iterrows(), start=1):
marker = " " if rank > 5 else "* "
print(f" {marker}{rank:<6} {name:<35} {row['mean']:.4f}")
print("\n* = included in plot")
print("="*60)
# %% [markdown]
# ### **t-SNE**
# %%
palette = 'tab10'
# %%
tSNE_X_mEPSC = mEPSC_even.drop("label", axis=1)
tSNE_y_mEPSC = mEPSC_even["label"]
# %%
tSNE_X_passive = passive_even.drop("label", axis=1)
tSNE_y_passive = passive_even["label"]
# %%
tSNE_X_spiking = spiking_even.drop("label", axis=1)
tSNE_y_spiking = spiking_even["label"]
# %%
# Re-run t-SNE with 3 components for mEPSC data
tSNE_mEPSC_3d = TSNE(n_components=3, verbose=1, perplexity=25, n_iter=300, random_state=42)
mEPSC_results_3d = tSNE_mEPSC_3d.fit_transform(tSNE_X_mEPSC)
mEPSC_tSNE_comp1_3d = mEPSC_results_3d[:,0]
mEPSC_tSNE_comp2_3d = mEPSC_results_3d[:,1]
mEPSC_tSNE_comp3_3d = mEPSC_results_3d[:,2]
# %%
# Re-run t-SNE with 3 components for passive data
tSNE_passive_3d = TSNE(n_components=3, verbose=1, perplexity=25, n_iter=300, random_state=42)
passive_results_3d = tSNE_passive_3d.fit_transform(tSNE_X_passive)
passive_tSNE_comp1_3d = passive_results_3d[:,0]
passive_tSNE_comp2_3d = passive_results_3d[:,1]
passive_tSNE_comp3_3d = passive_results_3d[:,2]
# %%
# Re-run t-SNE with 3 components for spiking data
tSNE_spiking_3d = TSNE(n_components=3, verbose=1, perplexity=25, n_iter=300, random_state=42)
spiking_results_3d = tSNE_spiking_3d.fit_transform(tSNE_X_spiking)
spiking_tSNE_comp1_3d = spiking_results_3d[:,0]
spiking_tSNE_comp2_3d = spiking_results_3d[:,1]
spiking_tSNE_comp3_3d = spiking_results_3d[:,2]
# %% [markdown]
# ### **Figures**
# %% [markdown]
# #### **Model vs Accuracy Figure**
# %%
#POSTER VERSION
plot_colors = {
"Action Potentials": "#4393c3",
"Passive": "#92c5de",
"mEPSC": "#2166ac"
}
model_labels = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
# Narrower width compresses points together; taller height gives more breathing room
fig, ax = plt.subplots(figsize=(10, 5), dpi=300)
plt.style.use('seaborn-v0_8-darkgrid')
# Lighter background
ax.set_facecolor('#f7f9fc')
fig.patch.set_facecolor('#f7f9fc')
datasets = {
"Action Potentials": (spiking_accuracies, spiking_std),
"Passive": (passive_accuracies, passive_std),
"mEPSC": (mEPSC_accuracies, mEPSC_std),
}
for label, (means, stds) in datasets.items():
ax.errorbar(
model_labels, means, yerr=stds,
label=label,
color=plot_colors[label],
marker='o',
linewidth=2.0,
markersize=6,
capsize=4,
capthick=1.5,
elinewidth=1.5,
)
ax.set_ylabel("Classification Accuracy", fontsize=16, fontweight='bold')
ax.set_xlabel("Model Type", fontsize=16, fontweight='bold')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.xticks(rotation=45, ha='right', fontsize=13, fontweight='bold')
plt.yticks(fontsize=13, fontweight='bold')
legend = ax.legend(
title="Feature Set",
fontsize=12,
title_fontsize=13,
frameon=False,
bbox_to_anchor=(1.02, 1.0),
loc='upper left'
)
# Bold legend title
legend.get_title().set_fontweight('bold')
plt.tight_layout()
plt.show()
# %%
#HELPER FUNCTIONS FOR PAPER VERSION
#Loading stats results for significance markers
tukey_df = pd.read_csv(base_path + 'MC for accuracies.csv', index_col=0)
tukey_df[['feat1','model1']] = tukey_df['group1'].str.split('_', n=1, expand=True)
tukey_df[['feat2','model2']] = tukey_df['group2'].str.split('_', n=1, expand=True)
model_order = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
def cld_within_feature(feature_set):
"""Compact letters comparing models WITHIN one feature set (for panels B/C/D)"""
sub = tukey_df[(tukey_df['feat1']==feature_set) & (tukey_df['feat2']==feature_set)]
G = nx.Graph()
G.add_nodes_from(model_order)
for _, row in sub.iterrows():
if row['reject']:
G.add_edge(row['model1'], row['model2'])
coloring = nx.coloring.greedy_color(G, strategy='largest_first')
colors_sorted = sorted(set(coloring.values()))
letter_map = {c: chr(97+i) for i, c in enumerate(colors_sorted)}
return {m: letter_map[coloring[m]] for m in model_order}
def cld_within_model(model):
"""Compact letters comparing feature sets WITHIN one model (for panel A)"""
sub = tukey_df[(tukey_df['model1']==model) & (tukey_df['model2']==model)]
G = nx.Graph()
G.add_nodes_from(['AP','Passive','mEPSC'])
for _, row in sub.iterrows():
if row['reject']:
G.add_edge(row['feat1'], row['feat2'])
coloring = nx.coloring.greedy_color(G, strategy='largest_first')
colors_sorted = sorted(set(coloring.values()))
letter_map = {c: chr(97+i) for i, c in enumerate(colors_sorted)}
return {f: letter_map[coloring[f]] for f in ['AP','Passive','mEPSC']}
def tukey_row(a, b):
"""Return the Tukey row for groups a,b in either order."""
r = tuk[((tuk.group1 == a) & (tuk.group2 == b)) |
((tuk.group1 == b) & (tuk.group2 == a))]
if len(r) != 1:
raise KeyError(f"expected 1 row for {a} vs {b}, got {len(r)}")
return r.iloc[0]
SIG_LEVELS = [(1e-4, '**'), (0.05, '*')]
def stars(p):
for thresh, sym in SIG_LEVELS:
if p < thresh:
return sym
return None
def _summarise(d):
mu = [np.mean(v) * 100 for v in d.values()]
sd = [np.std(v) * 100 for v in d.values()]
return mu, sd
def compact_letters(fs):
ns = {m: set() for m in model_labels}
for m1, m2 in itertools.combinations(model_labels, 2):
if stars(tukey_row(f'{fs}_{m1}', f'{fs}_{m2}')['p-adj']) is None:
ns[m1].add(m2)
ns[m2].add(m1)
order = sorted(model_labels, key=lambda m: -ACC[fs][model_labels.index(m)])
cliques = []
for m in order:
cand = [m] + [x for x in order if x in ns[m]]
for start in range(len(cand)):
cl = {cand[start]}
for x in cand:
if x in cl:
continue
if all(y in ns[x] for y in cl):
cl.add(x)
if cl not in cliques:
cliques.append(cl)
cliques = [c for c in cliques if not any(c < d for d in cliques)]
cliques.sort(key=lambda c: -max(ACC[fs][model_labels.index(m)] for m in c))
lab = {m: '' for m in model_labels}
for letter, c in zip('abcdefghij', cliques):
for m in c:
lab[m] += letter
return lab, cliques
def draw_panel(axis, fs):
a, e = ACC[fs], STD[fs]
axis.errorbar(xpos, a, yerr=e, color=plot_colors[SET2COL[fs]], marker='o',
linewidth=0, markersize=6, capsize=4, capthick=1.2, elinewidth=1.2)
for m1, m2 in itertools.combinations(model_labels, 2):
r = tukey_row(f'{fs}_{m1}', f'{fs}_{m2}')
audit.append((SET2COL[fs], f'{fs}: {m1} vs {m2}', int(r.csv_line), r.group1,
r.group2, r.meandiff, r['p-adj'], stars(r['p-adj']) or 'ns'))
top = 105
if ANNOT_MODE == 'cld':
lab, _ = compact_letters(fs)
for i, m in enumerate(model_labels):
axis.text(i, a[i] + e[i] + 2.0, ','.join(lab[m]), ha='center',
fontsize=11, fontweight='bold')
print(f"[{SET2COL[fs]}] CLD: " + ", ".join(
f"{m}={','.join(lab[m])}" for m in sorted(model_labels,
key=lambda m: -a[model_labels.index(m)])))
elif ANNOT_MODE == 'brackets':
sig = []
for m1, m2 in itertools.combinations(model_labels, 2):
sym = stars(tukey_row(f'{fs}_{m1}', f'{fs}_{m2}')['p-adj'])
if sym:
sig.append((model_labels.index(m1), model_labels.index(m2), sym))
sig.sort(key=lambda t: abs(t[1] - t[0]))
y = max(np.array(a) + np.array(e)) + 3
for i1, i2, sym in sig:
axis.plot([i1, i1, i2, i2], [y, y + 0.9, y + 0.9, y], color='black', lw=1)
axis.text((i1 + i2) / 2, y + 1.1, sym, ha='center', fontsize=7)
y += 3.4
top = y + 3
axis.set_ylim(40, top)
#CONSTANTS AND VARIABLES
letters_A = {m: cld_within_model(m) for m in model_order}
letters_AP = cld_within_feature('AP')
letters_Passive = cld_within_feature('Passive')
letters_mEPSC = cld_within_feature('mEPSC')
TUKEY_CSV = base_path + 'MC for accuracies.csv'
OUT_TIF = base_path + 'figure3.tif'
ANNOT_MODE = 'cld'
SHOW_NS_A = False
LABEL_ALL = True
DPI = 300
FIGSIZE = (16, 10)
model_labels = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
feature_sets = ['AP', 'Passive', 'mEPSC']
plot_colors = {"Action Potentials": "#4393c3", "Passive": "#92c5de", "mEPSC": "#2166ac"}
SET2COL = {'AP': "Action Potentials", 'Passive': "Passive", 'mEPSC': "mEPSC"}
ARM_START = 0.15
LANE_STEP = 0.26
ARM_LEN = 0.09
STAR_SIZE = 12
STAR_PAD = 0.03
STAR_STACK = True
# %%
tuk = pd.read_csv(TUKEY_CSV, index_col=0)
tuk['csv_line'] = tuk.index + 2 # 1-indexed line in the raw file
xpos = np.arange(len(model_labels))
XLIM = (-0.5, len(model_labels) + 0.05) # same on every panel so columns align
spiking_accuracies, spiking_std = _summarise(spiking_seed_accuracies)
passive_accuracies, passive_std = _summarise(passive_seed_accuracies)
mEPSC_accuracies, mEPSC_std = _summarise(mEPSC_seed_accuracies)
ACC = {'AP': spiking_accuracies, 'Passive': passive_accuracies, 'mEPSC': mEPSC_accuracies}
STD = {'AP': spiking_std, 'Passive': passive_std, 'mEPSC': mEPSC_std}
plt.style.use('seaborn-v0_8-darkgrid')
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=FIGSIZE)
audit = []
# ============ PANEL A : within model, across feature sets ============
for fs in feature_sets:
ax[0, 0].errorbar(xpos, ACC[fs], yerr=STD[fs], label=SET2COL[fs],
color=plot_colors[SET2COL[fs]], marker='o', linewidth=0,
markersize=5, capsize=4, capthick=1.2, elinewidth=1.2)
for i, m in enumerate(model_labels):
pairs = []
for a, b in itertools.combinations(feature_sets, 2):
r = tukey_row(f'{a}_{m}', f'{b}_{m}')
sym = stars(r['p-adj'])
audit.append(('A', f'{m}: {a} vs {b}', int(r.csv_line), r.group1, r.group2,
r.meandiff, r['p-adj'], sym or 'ns'))
if sym is None and not SHOW_NS_A:
continue
pairs.append((ACC[a][i], ACC[b][i], sym or 'ns'))
pairs.sort(key=lambda t: abs(t[0] - t[1]))
for lane, (y1, y2, sym) in enumerate(pairs):
x0 = i + ARM_START + lane * LANE_STEP
x1 = x0 + ARM_LEN
lo, hi = sorted([y1, y2])
ax[0, 0].plot([x1, x1], [lo, hi], color='black', lw=1)
ax[0, 0].plot([x0, x1], [lo, lo], color='black', lw=1)
ax[0, 0].plot([x0, x1], [hi, hi], color='black', lw=1)
txt = '\n'.join(sym) if (STAR_STACK and sym != 'ns') else sym
ax[0, 0].text(x1 + STAR_PAD, (lo + hi) / 2, txt, fontsize=STAR_SIZE,
style='italic' if sym == 'ns' else 'normal',
va='center', ha='left', linespacing=0.25)
ax[0, 0].set_ylim(40, 105)
ax[0, 0].legend(title="Feature Set", fontsize=9, title_fontsize=10,
frameon=False, loc='lower left')
draw_panel(ax[0, 1], 'AP')
draw_panel(ax[1, 0], 'Passive')
draw_panel(ax[1, 1], 'mEPSC')
# ============ shared cosmetics ============
uniform_y = len({tuple(a.get_ylim()) for a in ax.ravel()}) == 1
for (row, col), axis in np.ndenumerate(ax):
axis.set_xticks(xpos)
axis.set_xlim(*XLIM)
for sp in ('top', 'right'):
axis.spines[sp].set_visible(False)
show_y = LABEL_ALL or col == 0
show_x = LABEL_ALL or row == 1
axis.set_ylabel("Classification Accuracy (%)" if show_y else "", fontsize=12)
if not show_y and uniform_y:
axis.tick_params(labelleft=False)
axis.set_xlabel("Model Type" if show_x else "", fontsize=12)
if show_x:
axis.set_xticklabels(model_labels, rotation=45)
else:
axis.set_xticklabels([])
for axis, letter in zip(ax.ravel(), ['A', 'B', 'C', 'D']):
axis.text(-0.06, 1.06, letter, transform=axis.transAxes,
fontsize=15, fontweight='bold', va='top', ha='right')
plt.tight_layout()
fig.subplots_adjust(wspace=0.10, hspace=0.18)
# TIFF export. matplotlib renders RGBA; most journals require flattened RGB,
# so re-open and drop the alpha channel over white.
fig.savefig(OUT_TIF, dpi=DPI, facecolor='white',
pil_kwargs={'compression': 'tiff_lzw'})
_im = Image.open(OUT_TIF)
if _im.mode != 'RGB':
_flat = Image.new('RGB', _im.size, 'white')
_flat.paste(_im, mask=_im.split()[-1] if _im.mode == 'RGBA' else None)
_flat.save(OUT_TIF, compression='tiff_lzw', dpi=(DPI, DPI))
_im.close()
plt.show()
# ---------------- audit table for the supplement ----------------
audit_df = pd.DataFrame(audit, columns=['panel', 'comparison', 'csv_line',
'group1', 'group2', 'meandiff',
'p_adj', 'symbol'])
audit_df.to_csv(base_path + 'figure3_significance_audit.csv', index=False)
# %%
#Generated by Claude to assist with organizing MCs into output format for supplementary table S1
SUPP_PREFIX = base_path + 'figure3_supp_'
PP = 100.0
def fmt_p(p):
"""Adjusted p-values hit the solver's numerical floor (~1.9e-14), so anything
below the reporting threshold is shown as an inequality rather than a
spuriously precise number."""
if p < 1e-4:
return '<0.0001'
if p < 1e-3:
return f'{p:.5f}'
return f'{p:.4f}'
def oriented(g1, g2):
"""Tukey rows store an arbitrary group order. Return diff and CI expressed as
(mean of g2) - (mean of g1) regardless of how the row is stored."""
r = tukey_row(g1, g2)
diff, lo, hi = r.meandiff, r.lower, r.upper
if r.group1 != g1: # row is stored the other way round
diff, lo, hi = -diff, -hi, -lo # negate and swap the bounds
return diff * PP, lo * PP, hi * PP, r['p-adj'], int(r.csv_line)
# ---------------- Table S1: group summary ----------------
letters = {fs: compact_letters(fs)[0] for fs in feature_sets}
_seed_dicts = {'AP': spiking_seed_accuracies,
'Passive': passive_seed_accuracies,
'mEPSC': mEPSC_seed_accuracies}
rows_s1 = []
for fs in feature_sets:
for j, m in enumerate(model_labels):
rows_s1.append({
'Feature set': SET2COL[fs],
'Model': m,
'Mean accuracy (%)': round(ACC[fs][j], 2),
'SD (%)': round(STD[fs][j], 2),
'n (seeds)': len(list(_seed_dicts[fs].values())[j]),
'Group': ','.join(letters[fs][m]),
})
table_s1 = pd.DataFrame(rows_s1)
# ---------------- Table S2: pairwise comparisons ----------------
rows_s2 = []
# Panel A -- within model, across feature sets
for m in model_labels:
for a, b in itertools.combinations(feature_sets, 2):
d, lo, hi, p, line = oriented(f'{a}_{m}', f'{b}_{m}')
rows_s2.append({
'Panel': 'A', 'Stratum': m,
'Comparison': f'{SET2COL[b]} - {SET2COL[a]}',
'Difference (pp)': round(d, 2),
'95% CI': f'[{lo:.2f}, {hi:.2f}]',
'p-adj': fmt_p(p),
'Significance': stars(p) or 'n.s.',
'Source row': line,
})
# Panels B/C/D -- within feature set, across models
for fs, panel in zip(feature_sets, ['B', 'C', 'D']):
for m1, m2 in itertools.combinations(model_labels, 2):
d, lo, hi, p, line = oriented(f'{fs}_{m1}', f'{fs}_{m2}')
rows_s2.append({
'Panel': panel, 'Stratum': SET2COL[fs],
'Comparison': f'{m2} - {m1}',
'Difference (pp)': round(d, 2),
'95% CI': f'[{lo:.2f}, {hi:.2f}]',
'p-adj': fmt_p(p),
'Significance': stars(p) or 'n.s.',
'Source row': line,
})
table_s2 = pd.DataFrame(rows_s2)
# ---------------- sanity checks ----------------
assert len(table_s2) == 21 + 3 * 21, f"expected 84 comparisons, got {len(table_s2)}"
_recon = []
for _, r in table_s2.iterrows():
lo, hi = [float(v) for v in r['95% CI'].strip('[]').split(',')]
_recon.append(lo <= r['Difference (pp)'] <= hi)
assert all(_recon), "a difference fell outside its own CI -- check orientation logic"
table_s1.to_csv(SUPP_PREFIX + 'tableS1_group_summary.csv', index=False)
table_s2.to_csv(SUPP_PREFIX + 'tableS2_pairwise_tukey.csv', index=False)
# %%
print("="*10, "Model Accuracies and Standard Deviations", "="*10)
print("="*15, "mEPSC", "="*15)
for model, accuracies in mEPSC_seed_accuracies.items():
mu = np.mean(accuracies)
sigma = np.std(accuracies)
print(f"{model}: {mu:.4f} ± {sigma:.4f}")
print("="*15, "AP", "="*15)
for model, accuracies in spiking_seed_accuracies.items():
mu = np.mean(accuracies)
sigma = np.std(accuracies)
print(f"{model}: {mu:.4f} ± {sigma:.4f}")
print("="*15, "Passive", "="*15)
for model, accuracies in passive_seed_accuracies.items():
mu = np.mean(accuracies)
sigma = np.std(accuracies)
print(f"{model}: {mu:.4f} ± {sigma:.4f}")
# %% [markdown]
# #### **Classification Report Figure**
# %%
passive_diestrus_recalls = []
passive_estrus_recalls = []
passive_Eproestrus_recalls = []
passive_Lproestrus_recalls = []
for model_report in passive_classification_reports.values():
passive_diestrus_recalls.append(model_report['0']['recall'])
passive_estrus_recalls.append(model_report['1']['recall'])
passive_Eproestrus_recalls.append(model_report['2']['recall'])
passive_Lproestrus_recalls.append(model_report['3']['recall'])
spiking_diestrus_recalls = []
spiking_estrus_recalls = []
spiking_Eproestrus_recalls = []
spiking_Lproestrus_recalls = []
for model_report in spiking_classification_reports.values():
spiking_diestrus_recalls.append(model_report['0']['recall'])
spiking_estrus_recalls.append(model_report['1']['recall'])
spiking_Eproestrus_recalls.append(model_report['2']['recall'])
spiking_Lproestrus_recalls.append(model_report['3']['recall'])
mEPSC_diestrus_recalls = []
mEPSC_estrus_recalls = []
mEPSC_Eproestrus_recalls = []
mEPSC_Lproestrus_recalls = []
for model_report in mEPSC_classification_reports.values():
mEPSC_diestrus_recalls.append(model_report['0']['recall'])
mEPSC_estrus_recalls.append(model_report['1']['recall'])
mEPSC_Eproestrus_recalls.append(model_report['2']['recall'])
mEPSC_Lproestrus_recalls.append(model_report['3']['recall'])
# %%
print("Passive")
print(f"Diestrus: {np.mean(passive_diestrus_recalls):.4f} +- {np.std(passive_diestrus_recalls):.4f}")
print(f"Estrus: {np.mean(passive_estrus_recalls):.4f} +- {np.std(passive_estrus_recalls):.4f}")
print(f"Early Proestrus: {np.mean(passive_Eproestrus_recalls):.4f} +- {np.std(passive_Eproestrus_recalls):.4f}")
print(f"Late Proestrus: {np.mean(passive_Lproestrus_recalls):.4f} +- {np.std(passive_Lproestrus_recalls):.4f}")
print("AP")
print(f"Diestrus: {np.mean(spiking_diestrus_recalls):.4f} +- {np.std(spiking_diestrus_recalls):.4f}")
print(f"Estrus: {np.mean(spiking_estrus_recalls):.4f} +- {np.std(spiking_estrus_recalls):.4f}")
print(f"Early Proestrus: {np.mean(spiking_Eproestrus_recalls):.4f} +- {np.std(spiking_Eproestrus_recalls):.4f}")
print(f"Late Proestrus: {np.mean(spiking_Lproestrus_recalls):.4f} +- {np.std(spiking_Lproestrus_recalls):.4f}")
print("mEPSC")
print(f"Diestrus: {np.mean(mEPSC_diestrus_recalls):.4f} +- {np.std(mEPSC_diestrus_recalls):.4f}")
print(f"Estrus: {np.mean(mEPSC_estrus_recalls):.4f} +- {np.std(mEPSC_estrus_recalls):.4f}")
print(f"Early Proestrus: {np.mean(mEPSC_Eproestrus_recalls):.4f} +- {np.std(mEPSC_Eproestrus_recalls):.4f}")
print(f"Late Proestrus: {np.mean(mEPSC_Lproestrus_recalls):.4f} +- {np.std(mEPSC_Lproestrus_recalls):.4f}")
# %%
plt.style.use('seaborn-v0_8-white')
phase_labels = ['Diestrus', 'Estrus', 'Early Proestrus', 'Late Proestrus']
def make_long_df(diestrus, estrus, eproestrus, lproestrus, dataset_name):
model_names = list(passive_classification_reports.keys())
rows = []
for recalls, phase in zip([diestrus, estrus, eproestrus, lproestrus], phase_labels):
for model, val in zip(model_names, recalls):
rows.append({'phase': phase, 'Recall': val*100, 'model': model, 'Dataset': dataset_name})
return pd.DataFrame(rows)
passive_df = make_long_df(passive_diestrus_recalls, passive_estrus_recalls,
passive_Eproestrus_recalls, passive_Lproestrus_recalls, 'Passive')
spiking_df = make_long_df(spiking_diestrus_recalls, spiking_estrus_recalls,
spiking_Eproestrus_recalls, spiking_Lproestrus_recalls, 'AP')
mEPSC_df = make_long_df(mEPSC_diestrus_recalls, mEPSC_estrus_recalls,
mEPSC_Eproestrus_recalls, mEPSC_Lproestrus_recalls, 'mEPSC')
combined_df = pd.concat([spiking_df, passive_df, mEPSC_df], ignore_index=True)
class_colors = ['#E63946', '#457B9D', '#2A9D8F', '#E76F51']
fig, ax = plt.subplots(figsize=(8, 6), dpi=300)
fig.patch.set_facecolor('white')
ax.set_facecolor('white')
sns.stripplot(data=combined_df, x='Dataset', y='Recall', hue='phase',
hue_order=phase_labels, dodge=True, size=8, jitter=0.1,
alpha=0.8, palette=class_colors, ax=ax)
sns.pointplot(data=combined_df, x='Dataset', y='Recall', hue='phase',
hue_order=phase_labels, dodge=0.8 - 0.8/len(phase_labels),
estimator='mean', errorbar=None,
color='black', linestyle='none',
markers='_', markersize=20, markeredgewidth=3, ax=ax)
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles[:4], labels[:4], title='Phase', bbox_to_anchor=(1.05, 1),
loc='upper left', frameon=True, edgecolor='black', facecolor='white',
borderaxespad=0.)
sig_pairs = {
'AP': [((0, 1), '**'), ((2, 3), '***'), ((1, 3), '****')],
'Passive': [((0, 1), '*'), ((1, 3), '***')],
'mEPSC': [((2, 3), '*'), ((1, 3), '***'), ((0, 3), '*')],
}
dataset_order = ['AP', 'Passive', 'mEPSC']
n_hue = len(phase_labels)
hue_offset = {i: -0.4 + (0.8 / n_hue) * (i + 0.5) for i in range(n_hue)} # -0.3,-0.1,0.1,0.3
tick_h = 2.5
pad = 6.0
step = 9.0
for ds, pairs in sig_pairs.items():
xc = dataset_order.index(ds)
y_base = combined_df.loc[combined_df['Dataset'] == ds, 'Recall'].max() + pad
for k, ((a, b), stars) in enumerate(pairs):
y = y_base + k * step
x1 = xc + hue_offset[a]
x2 = xc + hue_offset[b]
ax.plot([x1, x1, x2, x2], [y, y + tick_h, y + tick_h, y],
lw=1.5, c='black', solid_capstyle='projecting', clip_on=False)
ax.text((x1 + x2) / 2, y + tick_h, stars, ha='center', va='bottom',
color='black', fontsize=13)
ax.set_ylim(0, 130)
ax.set_ylabel("Recall (%)", fontsize=14)
ax.set_xlabel("Dataset", fontsize=14)
fig.subplots_adjust(right=0.75)
plt.tight_layout()
# Explicitly set y-ticks from 0 to 100 even though the graph extends above that for significance bars
custom_yticks = np.arange(0, 101, 20)
ax.set_yticks(custom_yticks)
ax.set_yticklabels([str(int(t)) for t in custom_yticks])
plt.savefig(base_path + 'figure6.tif', bbox_inches='tight')
plt.show()
# %% [markdown]
# #### **Feature Importances Figure**
# %%
# Set the plot style and color
plt.style.use('seaborn-v0_8-darkgrid')
plot_color = "#2166ac"
# Create figure with 3 subplots
fig, axes = plt.subplots(3, 1, figsize=(4, 8), dpi=300)
# ============== Panel A: Spiking Dataset ==============
# Top 5 features
axes[0].bar(
range(5),
spiking_importances["mean"]*100,
color=plot_color,
linewidth=0
)
axes[0].set_ylabel("Feature importance (%)", fontsize=7)
axes[0].set_title("A", fontsize=10, fontweight='bold', loc='left', x=-0.125, y=1)
axes[0].set_xticks(range(5))
axes[0].set_xticklabels(spiking_top5, rotation=45, ha='right', fontsize=7)
axes[0].tick_params(axis='y', labelsize=7)
axes[0].spines['top'].set_visible(False)
axes[0].spines['right'].set_visible(False)
# ============== Panel B: Passive Dataset ==============
# Top 5 features
axes[1].bar(
range(5),
passive_importances["mean"]*100,
color=plot_color,
linewidth=0
)
axes[1].set_ylabel("Feature importance (%)", fontsize=7)
axes[1].set_title("B", fontsize=10, fontweight='bold', loc='left', x=-0.125, y=1)
axes[1].set_xticks(range(5))
axes[1].set_xticklabels(passive_top5, rotation=45, ha='right', fontsize=7)
axes[1].tick_params(axis='y', labelsize=7)
axes[1].spines['top'].set_visible(False)
axes[1].spines['right'].set_visible(False)
# ============== Panel C: mEPSC Dataset ==============
# Top 5 features
axes[2].bar(
range(5),
mEPSC_importances["mean"]*100,
color=plot_color,
linewidth=0
)
axes[2].set_ylabel("Feature importance (%)", fontsize=7)
axes[2].set_title("C", fontsize=10, fontweight='bold', loc='left', x=-0.125, y=1)
axes[2].set_xticks(range(5))
axes[2].set_xticklabels(mEPSC_top5, rotation=45, ha='right', fontsize=7)
axes[2].tick_params(axis='y', labelsize=7)
axes[2].spines['top'].set_visible(False)
axes[2].spines['right'].set_visible(False)
# Set the same y-axis limits for all three panels
for ax in axes:
ax.set_ylim(0, 35) # Add 5% padding at the top
plt.tight_layout()
plt.savefig(base_path + 'figure4.svg', bbox_inches='tight')
plt.show()
# %% [markdown]
# #### **Data Types Figure**
# %%
fig2_X_mEPSC_diestrus = read_abf(raw_mEPSC_diestrus[0])
fig2_X_mEPSC_estrus = read_abf(raw_mEPSC_estrus[0])
fig2_X_mEPSC_Eproestrus = read_abf(raw_mEPSC_Eproestrus[0])
fig2_X_mEPSC_Lproestrus = read_abf(raw_mEPSC_Lproestrus[0])
fig2_X_spike_diestrus = read_abf(raw_spike_diestrus[0])
fig2_X_spike_estrus = read_abf(raw_spike_estrus[0])
fig2_X_spike_Eproestrus = read_abf(raw_spike_Eproestrus[0])
fig2_X_spike_Lproestrus = read_abf(raw_spike_Lproestrus[0])
fig2_X_passive_diestrus = read_abf(raw_passive_diestrus[0])
fig2_X_passive_estrus = read_abf(raw_passive_estrus[0])
fig2_X_passive_Eproestrus = read_abf(raw_passive_Eproestrus[0])
fig2_X_passive_Lproestrus = read_abf(raw_passive_Lproestrus[0])
# %%
abf3 = pyabf.ABF(raw_spike_diestrus[0])
abf3.dataPointsPerMs*1000
# %%
abf = pyabf.ABF(raw_mEPSC_diestrus[0])
abf.dataPointsPerMs*1000
# %%
def filter_signal(current, freq=20000):
raw_signal = current
adjusted_signal = raw_signal - np.median(raw_signal)
b_lowpass, a_lowpass = signal.bessel(4, 800, 'low', analog=False, norm='phase', fs=freq)
b_notch, a_notch = signal.iirnotch(60.0, 30.0, fs=freq)
b_multiband = signal.convolve(b_lowpass, b_notch)
a_multiband = signal.convolve(a_lowpass, a_notch)
filtered_signal = signal.filtfilt(b_multiband, a_multiband, adjusted_signal)
return filtered_signal
# %%
filtered = filter_signal(fig2_X_mEPSC_diestrus)
# %%
abf4 = pyabf.ABF(raw_spike_diestrus[0])
x = len(fig2_X_mEPSC_diestrus)/abf4.sweepCount
# %%
sig = fig2_X_spike_diestrus[800000:835000]
# %%
abf2 = pyabf.ABF(raw_passive_diestrus[0])
# %%
len(abf2.sweepX)
# %%
# Data figure - 3 side by side plots
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(16, 12), dpi=300) # Changed to 2x2 grid and adjusted figsize
plt.style.use('seaborn-v0_8-darkgrid')
plot_color = "#2166ac"
# Panel A: Spiking plot
spike_data = fig2_X_spike_diestrus[800000:835000]
time_spike = np.arange(len(spike_data)) / abf3.dataPointsPerMs
ax[0,0].plot(time_spike, spike_data, color=plot_color, linewidth=1.0)
ax[0,0].set_ylabel("Voltage (mV)", fontsize=13)
ax[0,0].set_xlabel("Time (ms)", fontsize=13)
ax[0,0].spines['top'].set_visible(False)
ax[0,0].spines['right'].set_visible(False)
# Panel B: Sweep data
for sweep in abf2.sweepList:
abf2.setSweep(sweep)
time_sweep = abf2.sweepX[:40000] * 1000
ax[0,1].plot(time_sweep, abf2.sweepY[:40000], color=plot_color,
linewidth=0.5, alpha=0.7)
ax[0,1].set_ylabel("Voltage (mV)", fontsize=13)
ax[0,1].set_xlabel("Time (ms)", fontsize=13)
ax[0,1].spines['top'].set_visible(False)
ax[0,1].spines['right'].set_visible(False)
# Panel C: Filtered current
current_data = filtered[235000:255000]
time_current = np.arange(len(current_data)) / abf.dataPointsPerMs
ax[1,0].plot(time_current, current_data, color=plot_color, linewidth=1.0)
ax[1,0].set_ylabel("Current (pA)", fontsize=13)
ax[1,0].set_xlabel("Time (ms)", fontsize=13)
ax[1,0].spines['top'].set_visible(False)
ax[1,0].spines['right'].set_visible(False)
#Panel D: Burst metrics
last_burst = spiking_diestrus.iloc[1]
spike_data = fig2_X_spike_diestrus[800000:835000]
time_spike = np.arange(len(spike_data)) / abf3.dataPointsPerMs
ax[1,1].plot(time_spike, spike_data, color=plot_color, linewidth=1.0)
ax[1,1].set_ylabel("Voltage (mV)", fontsize=13)
ax[1,1].set_xlabel("Time (ms)", fontsize=13)
ax[1,1].spines['top'].set_visible(False)
ax[1,1].spines['right'].set_visible(False)
# --- Panel D annotations ---
burst_start_ms = last_burst['burst_start'] * 1000 - (800000 / abf3.dataPointsPerMs)
burst_end_ms = last_burst['burst_end'] * 1000 - (800000 / abf3.dataPointsPerMs)
peak_val = last_burst['avg_spike_peaks'] # 29.38 mV
amp_val = last_burst['avg_spike_amplitude'] # 81.51 mV
base_val = peak_val - amp_val # true amplitude foot (~ -52 mV)
# --- 1. Burst duration ---
plateau_val = np.percentile(spike_data, 65)
y_duration = plateau_val
ax[1,1].annotate('', xy=(burst_end_ms + 30, y_duration), xytext=(burst_start_ms - 20, y_duration),
arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))
ax[1,1].text((burst_start_ms + burst_end_ms) / 2, y_duration - 3,
f"Burst duration\n{last_burst['burst_length']:.2f} s",
ha='center', va='top', fontsize=12)
# --- 2. Avg spike amplitude (bracket now spans peak -> peak - amplitude) ---
bracket_x = burst_start_ms - 200 # sits left of the rise, over flat baseline
tick_width = 8
ax[1,1].plot([bracket_x, bracket_x], [base_val + 5, peak_val], color='black', lw=1.2)
ax[1,1].plot([bracket_x, bracket_x + tick_width], [peak_val, peak_val], color='black', lw=1.2)
ax[1,1].plot([bracket_x, bracket_x + tick_width], [base_val + 5, base_val + 5], color='black', lw=1.2)
ax[1,1].text(bracket_x - 12, (peak_val + base_val) / 2,
f"Average spike\namplitude\n{amp_val:.2f} mV",
color='black', fontsize=12, va='center', ha='right')
# --- 3. Avg spike peak value (dashed line only across the spike train) ---
pad = 25
ax[1,1].plot([burst_start_ms - pad, burst_end_ms + pad], [peak_val, peak_val],
color='black', linestyle='--', linewidth=0.8)
ax[1,1].text(burst_end_ms + pad + 20, peak_val,
f"Average spike\npeak value\n{peak_val:.2f} mV",
ha='left', va='center', fontsize=12, color='black')
# --- 4. Spikes per burst ---
ax[1,1].text(burst_end_ms + pad + 20, -5,
f"Spikes per burst = 5",
ha='left', va='center', fontsize=12, color='black')
# Add panel labels
for i, label in enumerate(['A', 'B', 'C', 'D']):
row = i // 2
col = i % 2
ax[row, col].text(-0.05, 1.05, label, transform=ax[row, col].transAxes,
fontsize=16, fontweight='bold', va='top', ha='right')
plt.tight_layout()
plt.savefig(base_path + 'figure1.tif', bbox_inches='tight')
plt.show()
# %% [markdown]
# #### **t-SNE Figure**
# %%
plt.style.use('default')
views = [(30, 45), (15, 180), (45, 270)]
zoom = 1.2
data_sets = [
(spiking_tSNE_comp1_3d, spiking_tSNE_comp2_3d, spiking_tSNE_comp3_3d, tSNE_y_spiking),
(passive_tSNE_comp1_3d, passive_tSNE_comp2_3d, passive_tSNE_comp3_3d, tSNE_y_passive),
(mEPSC_tSNE_comp1_3d, mEPSC_tSNE_comp2_3d, mEPSC_tSNE_comp3_3d, tSNE_y_mEPSC)
]
row_titles = ['A', 'B', 'C']
label_style = {
(30, 45): {'x': {'pad': -6, 'size': 9}, 'y': {'pad': -6, 'size': 9}, 'z': {'pad': -10, 'size': 9}},
(15, 180): {'x': {'pad': -8, 'size': 7}, 'y': {'pad': -2, 'size': 10}, 'z': {'pad': -10, 'size': 9}},
(45, 270): {'x': {'pad': -6, 'size': 9}, 'y': {'pad': -6, 'size': 9}, 'z': {'pad': -6, 'size': 9}},
}
fig = plt.figure(figsize=(22, 18), dpi=300)
fig.patch.set_facecolor('white')
plt.subplots_adjust(wspace=0.15, hspace=0.15)
scatter = None
class_colors = ['#E63946', '#457B9D', '#2A9D8F', '#E76F51']
cmap = mcolors.ListedColormap(class_colors)
for row, (comp1, comp2, comp3, labels) in enumerate(data_sets):
for col, (elev, azim) in enumerate(views):
ax = fig.add_subplot(3, 3, row * 3 + col + 1, projection='3d')
ax.set_facecolor('white')
ax.set_box_aspect(None, zoom=zoom)
scatter = ax.scatter(
comp1, comp2, comp3,
c=labels,
cmap=cmap,
alpha=1,
s=6,
vmin=0, vmax=3,
edgecolor='w',
linewidth=0.25
)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_zticklabels([])
ax.tick_params(axis='both', which='both', length=0)
ax.view_init(elev=elev, azim=azim)
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False
ax.xaxis.pane.set_edgecolor('lightgrey')
ax.yaxis.pane.set_edgecolor('lightgrey')
ax.zaxis.pane.set_edgecolor('lightgrey')
ax.grid(True, color='lightgrey', linewidth=0.5)
if col == 0:
ax.set_title(row_titles[row], fontsize=20, fontweight='bold', loc='left', pad=2)
legend_labels = ['Diestrus', 'Estrus', 'Early proestrus', 'Late proestrus']
legend_handles = [
plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=c,
markersize=9, label=l)
for c, l in zip(class_colors, legend_labels)
]
fig.legend(
handles=legend_handles,
labels=legend_labels,
fontsize=20,
title_fontsize=0,
frameon=False,
loc='center right',
bbox_to_anchor=(1.02, 0.5),
handletextpad=0.5,
labelspacing=0.8,
)
plt.savefig(base_path + 'figure5.svg', bbox_inches='tight')
plt.show()
# %% [markdown]
# ### **Stats**
# %% [markdown]
# #### **AP vs. passive vs. mEPSC accuracies**
# %%
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
all_accuracies_data = []
# Populate with mEPSC seed accuracies (individual values from each seed)
for name in model_names:
for acc in mEPSC_seed_accuracies[name]:
all_accuracies_data.append({
"Accuracy": acc,
"Feature_Set": "mEPSC",
"Model": name
})
# Populate with spiking seed accuracies
for name in model_names:
for acc in spiking_seed_accuracies[name]:
all_accuracies_data.append({
"Accuracy": acc,
"Feature_Set": "AP",
"Model": name
})
# Populate with passive seed accuracies
for name in model_names:
for acc in passive_seed_accuracies[name]:
all_accuracies_data.append({
"Accuracy": acc,
"Feature_Set": "Passive",
"Model": name
})
accuracies_df = pd.DataFrame(all_accuracies_data)
# Drop rows with NaN or inf values in the 'Accuracy' column just in case
accuracies_df = accuracies_df.replace([np.inf, -np.inf], np.nan).dropna(subset=['Accuracy'])
formula = "Accuracy ~ C(Feature_Set) * C(Model)"
model = smf.ols(formula, data=accuracies_df).fit()
table = sm.stats.anova_lm(model, typ=2)
display(table)
# %%
df_grouped = accuracies_df.copy()
df_grouped['group'] = df_grouped['Feature_Set'] + '_' + df_grouped['Model']
tukey = pairwise_tukeyhsd(df_grouped['Accuracy'], df_grouped['group'])
tukey_df = pd.DataFrame(data=tukey._results_table.data[1:],
columns=tukey._results_table.data[0])
tukey_df['p-adj'] = tukey.pvalues #To get precise values
pd.set_option('display.float_format', lambda x: f'{x:.12f}')
tukey_df.to_csv("/content/drive/MyDrive/Estrous Classification Project/MC for accuracies.csv")
# %% [markdown]
# #### **Across full vs. subsetted dataset**
# %%
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
SEEDS = [42, 7, 13, 21, 99]
# Prepare data for ANOVA
all_accuracies_data = []
# Full dataset accuracies
for model_name in model_names:
for seed_acc in mEPSC_seed_accuracies[model_name]:
all_accuracies_data.append({
'Accuracy': seed_acc,
'Dataset_Type': 'Full',
'Feature_Set': 'mEPSC',
'Model': model_name
})
for seed_acc in spiking_seed_accuracies[model_name]:
all_accuracies_data.append({
'Accuracy': seed_acc,
'Dataset_Type': 'Full',
'Feature_Set': 'Spiking',
'Model': model_name
})
for seed_acc in passive_seed_accuracies[model_name]:
all_accuracies_data.append({
'Accuracy': seed_acc,
'Dataset_Type': 'Full',
'Feature_Set': 'Passive',
'Model': model_name
})
# Even label count dataset accuracies
for model_name in model_names:
for seed_acc in mEPSC_even_seed_accuracies[model_name]:
all_accuracies_data.append({
'Accuracy': seed_acc,
'Dataset_Type': 'Even_Label',
'Feature_Set': 'mEPSC',
'Model': model_name
})
for seed_acc in spiking_even_seed_accuracies[model_name]:
all_accuracies_data.append({
'Accuracy': seed_acc,
'Dataset_Type': 'Even_Label',
'Feature_Set': 'Spiking',
'Model': model_name
})
for seed_acc in passive_even_seed_accuracies[model_name]:
all_accuracies_data.append({
'Accuracy': seed_acc,
'Dataset_Type': 'Even_Label',
'Feature_Set': 'Passive',
'Model': model_name
})
accuracies_df = pd.DataFrame(all_accuracies_data)
formula = 'Accuracy ~ C(Dataset_Type) * C(Feature_Set) * C(Model)'
model = smf.ols(formula, data=accuracies_df).fit()
anova_table = sm.stats.anova_lm(model, typ=2) # Type 2 ANOVA for unbalanced designs
print("\nANOVA Results:")
display(anova_table)
# %%
df_grouped = accuracies_df.copy()
df_grouped['group'] = df_grouped['Dataset_Type'] + '_' + df_grouped['Feature_Set']
tukey = pairwise_tukeyhsd(df_grouped['Accuracy'], df_grouped['group'])
print(tukey.summary())
# %% [markdown]
# #### **Across randomized vs. non-randomized**
# %%
model_names = ['RFC', 'GBC', 'LR', 'KNN', 'MLP', 'SVC', 'DTC']
SEEDS = [42, 7, 13, 21, 99]
# Prepare data for ANOVA
all_accuracies_rnd_data = []
# Full dataset accuracies — real labels
for model_name in model_names:
for seed_acc in mEPSC_seed_accuracies[model_name]:
all_accuracies_rnd_data.append({
'Accuracy': seed_acc,
'Label_Type': 'Real',
'Feature_Set': 'mEPSC',
'Model': model_name
})
for seed_acc in spiking_seed_accuracies[model_name]:
all_accuracies_rnd_data.append({
'Accuracy': seed_acc,
'Label_Type': 'Real',
'Feature_Set': 'Spiking',
'Model': model_name
})
for seed_acc in passive_seed_accuracies[model_name]:
all_accuracies_rnd_data.append({
'Accuracy': seed_acc,
'Label_Type': 'Real',
'Feature_Set': 'Passive',
'Model': model_name
})
# Randomized label control accuracies
for model_name in model_names:
for seed_acc in mEPSC_rnd_seed_accuracies[model_name]:
all_accuracies_rnd_data.append({
'Accuracy': seed_acc,
'Label_Type': 'Randomized',
'Feature_Set': 'mEPSC',
'Model': model_name
})
for seed_acc in spiking_rnd_seed_accuracies[model_name]:
all_accuracies_rnd_data.append({
'Accuracy': seed_acc,
'Label_Type': 'Randomized',
'Feature_Set': 'Spiking',
'Model': model_name
})
for seed_acc in passive_rnd_seed_accuracies[model_name]:
all_accuracies_rnd_data.append({
'Accuracy': seed_acc,
'Label_Type': 'Randomized',
'Feature_Set': 'Passive',
'Model': model_name
})
accuracies_rnd_df = pd.DataFrame(all_accuracies_rnd_data)
formula = 'Accuracy ~ C(Label_Type) * C(Feature_Set) * C(Model)'
rnd_anova_model = smf.ols(formula, data=accuracies_rnd_df).fit()
anova_rnd_table = sm.stats.anova_lm(rnd_anova_model, typ=2) # Type 2 for unbalanced designs
display(anova_rnd_table)
# %% [markdown]
# #### **Across Phase and Feature Set in Classification Reports**
# %%
all_accuracies_data = []
for name in model_names:
rep = mEPSC_classification_reports[name]
iter = 0
for key, value in rep.items():
if iter > 3:
break
all_accuracies_data.append({
"Recall": value['recall'], #recall
"Model": name,
"Phase": key,
"Feature_Set": "mEPSC"
})
iter += 1
for name in model_names:
rep = spiking_classification_reports[name]
iter = 0
for key, value in rep.items():
if iter > 3:
break
all_accuracies_data.append({
"Recall": value['recall'], #recall
"Model": name,
"Phase": key,
"Feature_Set": "AP"
})
iter += 1
for name in model_names:
rep = passive_classification_reports[name]
iter = 0
for key, value in rep.items():
if iter > 3:
break
all_accuracies_data.append({
"Recall": value['recall'], #recall
"Model": name,
"Phase": key,
"Feature_Set": "Passive"
})
iter += 1
accuracies_df = pd.DataFrame(all_accuracies_data)
accuracies_df = accuracies_df.replace([np.inf, -np.inf], np.nan).dropna(subset=['Recall'])
formula = "Recall ~ C(Phase) * C(Feature_Set) + C(Phase) * C(Model) + C(Feature_Set) * C(Model)"
model = smf.ols(formula, data=accuracies_df).fit()
table = sm.stats.anova_lm(model, typ=2)
display(table)
# %%
df_grouped = accuracies_df.copy()
df_grouped['group'] = df_grouped['Phase'] + '_' + df_grouped['Feature_Set']
tukey = pairwise_tukeyhsd(df_grouped['Recall'], df_grouped['group'])
tukey_df = pd.DataFrame(data=tukey._results_table.data[1:],
columns=tukey._results_table.data[0])
tukey_df['p-adj'] = tukey.pvalues
pd.set_option('display.float_format', lambda x: f'{x:.12f}')
tukey_df.to_csv("/content/drive/MyDrive/Estrous Classification Project/MC for recall.csv")
print(tukey.summary())
# %%
df_grouped = accuracies_df.copy()
df_grouped['group'] = df_grouped['Feature_Set'] + '_' + df_grouped['Model']
tukey = pairwise_tukeyhsd(df_grouped['Accuracy'], df_grouped['group'])
tukey_df = pd.DataFrame(data=tukey._results_table.data[1:],
columns=tukey._results_table.data[0])
tukey_df['p-adj'] = tukey.pvalues
pd.set_option('display.float_format', lambda x: f'{x:.12f}')
tukey_df.to_csv("/content/drive/MyDrive/Estrous Classification Project/MC for accuracies.csv")

Estrous_Classification.ipynb at commit 8495b02, under MIT · at the source

Overview

Authors: Armaan Raina¹, John Meitzen^1,2,3

ORCID iDs: John Meitzen

Department of Biological Sciences, North Carolina State University, Raleigh, NC USA
Center for Human Health and the Environment, North Carolina State University, Raleigh, NC USA
Dept. of Biological Sciences, NC State University, 144 David Clark Labs, Campus Box 7617, Raleigh, NC 27695-7617 USA

Institutions: North Carolina State University (United States)

Journal: Neuroinformatics, volume 24, issue 3, article 57

Dates: received 8 July 2026; accepted 25 August 2026; published online 1 September 2026; in print 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1007/s12021-026-09814-0 · PMID 42678460 · PMCID PMC13534189 · OpenAlex W7204910555

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: intracellular / patch clamp (modality), rat (organism), methods / tools (subfield)

Methods: Smoothing, state filtering, decompositions, Statistics, Machine learning, Preprocessing, Evoked potentials, Single-unit activity, calcium imaging

Keywords: machine learning, estrous cycle, nucleus accumbens, electrophysiology, neuroinformatics

MeSH: Estrous Cycle*, Machine Learning*, Models, Neurological*, Neurons*, Nucleus Accumbens*, Action Potentials, Animals, Classification Algorithms, Excitatory Postsynaptic Potentials, Female, Medium Spiny Neurons, Predictive Learning Models, Rats (* major topic)

Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: NIEHS NIH HHS (P30ES025128, P30 ES025128); North Carolina State University (Provost Professional Experience Award); National Institute of Environmental Health Sciences (P30ES025128)

Citations: not cited yet (Europe PMC); 73 references in the paper

Abstract

Machine learning models (MLMs) have been used to classify neuron subtypes based upon electrophysiological attributes. What has been less explored is using MLMs to assess differences within a neuron type, especially in the context of subtle changes induced by neuromodulation such as the rodent estrous cycle. Previous research found that the estrous cycle shifts rat nucleus accumbens (NAc) medium spiny neuron (MSN) electrophysiology, including action potential (AP), passive, and miniature excitatory post-synaptic current (mEPSC) properties. This plasticity provides a model system to investigate this question. We hypothesized that MLM classification accuracy would differ when trained across these data types, reflecting the information each encodes regarding estrous phase identity. To test this hypothesis, we extracted electrophysiological features across four estrous cycle phases from a publicly available dataset, and employed this data to train MLMs to classify estrous cycle phase origin. We found: MLMs identified estrous phase origin with up to 94% accuracy (Random Forest); MLM accuracy differed by model and feature set; passive, AP and mEPSC features yielded sequential performances with feature importance analysis revealing input resistance, AP burst length, and mEPSC baseline current as most discriminative for their respective feature sets; the late proestrus phase demonstrated the most distinct electrophysiological profile. These findings demonstrate that MLMs can detect changes associated with the estrous cycle and are useful for assessing electrophysiological differences within a neuron type. They further represent a broader application of MLMs as computational tools useful for understanding the sensitivity of neural properties to the influences of neuromodulatory cycles.

Supplementary Information: The online version contains supplementary material available at https://doi.org/10.1007/s12021-026-09814-0.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

doi:10.5061/dryad.k0p2ngfgm

License: none: the authors keep all their rights

State: the link answers, verified on 26 September 2026

Evidence: the link answers

Software Heritage: not checked

Found in: “Data Availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)

26 September 2026: the link answers (HTTP 200)

At the source: datadryad.org/dataset/doi:10.5061/dryad.k0p2ngfgm

Armaan-Raina/Estrous-Phase-Classification

License: MIT

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Commit: 8495b02c0cd77faeb25c2218dff6c0b2a9311508, 25 August 2026

Languages: Jupyter (1)

Size: 3 files, 1 script

Software Heritage: not archived

Found in: “Data Availability”

Holds: README, license file, 1 notebook

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Matplotlib (1 file), NetworkX (1 file), NumPy (1 file), pandas (1 file), Pillow (1 file), pyABF (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)

Availability: 1 check, the latest on 26 September 2026: the link answers

26 September 2026: the link answers

3 files

Estrous_Classification.ipynb, Jupyter, 2,588 lines, 4 matches
LICENSE, License, 21 lines
README.md, Text, 15 lines

At the source: github.com/Armaan-Raina/Estrous-Phase-Classification

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

The parent electrophysiological dataset used in this study is available at the Dryad Data Repository (https://datadryad.org/dataset/doi:10.5061/dryad.k0p2ngfgm). The code generated during this study to analyze the parent dataset is available in the GitHub repository (https://github.com/Armaan-Raina/Estrous-Phase-Classification).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 3, 28 September 2026

Publisher: n/a → Springer Science+Business Media

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 2 authors, 5 keywords, 13 MeSH terms, 3 funders, 69 references.

Cite

This paper

Raina, A., & Meitzen, J. (2026). Application of Machine Learning Models to Identify Differences in Neural Electrophysiological Properties Across Estrous Cycle Phases. Neuroinformatics, 24(3), 57. https://doi.org/10.1007/s12021-026-09814-0

BibTeX

@article{raina2026application, author = {Raina, Armaan and Meitzen, John}, title = {{Application of Machine Learning Models to Identify Differences in Neural Electrophysiological Properties Across Estrous Cycle Phases}}, journal = {Neuroinformatics}, year = {2026}, month = sep, volume = {24}, number = {3}, pages = {57}, publisher = {Springer Science+Business Media}, issn = {1539-2791}, doi = {10.1007/s12021-026-09814-0}, url = {https://doi.org/10.1007/s12021-026-09814-0}, pmid = {42678460}, pmcid = {PMC13534189} }

RIS

TY - JOUR AU - Raina, Armaan AU - Meitzen, John TI - Application of Machine Learning Models to Identify Differences in Neural Electrophysiological Properties Across Estrous Cycle Phases T2 - Neuroinformatics J2 - Neuroinformatics PY - 2026 DA - 2026/09/01 VL - 24 IS - 3 SP - 57 SN - 1539-2791 PB - Springer Science+Business Media DO - 10.1007/s12021-026-09814-0 UR - https://doi.org/10.1007/s12021-026-09814-0 LA - en ER -

CSL-JSON

{ "id": "10.1007/s12021-026-09814-0", "type": "article-journal", "title": "Application of Machine Learning Models to Identify Differences in Neural Electrophysiological Properties Across Estrous Cycle Phases", "container-title": "Neuroinformatics", "author": [ { "family": "Raina", "given": "Armaan" }, { "family": "Meitzen", "given": "John" } ], "container-title-short": "Neuroinformatics", "volume": "24", "issue": "3", "page": "57", "DOI": "10.1007/s12021-026-09814-0", "PMID": "42678460", "PMCID": "PMC13534189", "ISSN": "1539-2791", "publisher": "Springer Science+Business Media", "URL": "https://doi.org/10.1007/s12021-026-09814-0", "language": "en", "issued": { "date-parts": [ [ 2026, 9, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-76939-w [code]: HIPPIE: a generative model for electrophysiological analysis across species, technologies, and modalities.
Journal: Nature communications
In common: Pillow, statsmodels, seaborn, 5 other tools, rat, methods / tools, 1 reference
[2] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: NetworkX, Pillow, statsmodels, 6 other tools
[3] doi:10.3390/ijms27167275 [code]: Integrative Multi-Omics Analysis of Multiple Sclerosis Reveals Cell-Type-Specific Regulatory Landscapes and Discordant Methylation-Expression Coupling.
Journal: International journal of molecular sciences
In common: NetworkX, Pillow, statsmodels, 6 other tools
[4] doi:10.1371/journal.pcbi.1014571 [code]: SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: NetworkX, Pillow, statsmodels, 6 other tools
[5] doi:10.1016/j.isci.2026.116825 [code]: Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: NetworkX, Pillow, statsmodels, 6 other tools
[6] doi:10.1038/s41467-026-74357-6 [code]: Hippocampo-neocortical interaction as compressive retrieval-augmented generation.
Journal: Nature communications
In common: NetworkX, Pillow, statsmodels, 6 other tools
[7] doi:10.1016/j.isci.2026.116206 [code]: Gut distension evokes rapid neural dynamics in vagal and hindbrain populations of larval zebrafish.
Journal: iScience
In common: NetworkX, Pillow, statsmodels, 6 other tools
[8] doi:10.3389/fncom.2026.1786996 [code]: Schumann-anchored golden ratio organization of human neural oscillations.
Journal: Frontiers in computational neuroscience
In common: NetworkX, Pillow, statsmodels, 6 other tools
[9] doi:10.1016/j.celrep.2026.117420 [code]: Neural population dynamics of direct electrical stimulation of neocortex.
Journal: Cell reports
In common: NetworkX, Pillow, statsmodels, 6 other tools
[10] doi: [code]: Naturalistic behavior and self-generated neural activity predictive of self-correction
Journal: bioRxiv : the preprint server for biology
In common: NetworkX, Pillow, statsmodels, 6 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.