OSCR

The role of amygdala GABA neurons in controlling stress and reproduction in female mice.

Code ↔ Paper

17 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 17 matches · 4 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Comparing pre-intervention and intervention cell-cell interaction strength ↔ data_analysis/calculate_riemannian_distances.ipynb, lines 12–127 · score 0.83 · Metropolis Hasting, SPD matrices, Riemannian distance, correlation matrices, symmetric, uniform
  2. [2] § Methods › Comparing pre-intervention and intervention cell-cell interaction strength ↔ data_analysis/analysis_pipeline.ipynb, lines 303–418 · score 0.82 · Metropolis Hasting, SPD matrices, Riemannian distance, correlation matrices, symmetric, uniform
  3. [3] § Results › In silico modeling reveals complementary relationship between MePD GABA and UCN3 populations ↔ MePD_model/simulations/KNDyXMePDU.m, lines 2–61 · score 0.71 · neuronal activity, KNDy model, model parameter, UCN3 stimulation, couple, secretion
  4. [4] § Results › Dynamic reconfiguration of MePD GABA neuronal networks during selective optogenetic stimulation of urocortin-3 (UCN3) neurons and restraint stress ↔ data_analysis/data/calcium_imaging/figures/norm_Riemannian distance_stat_fig_1.ipynb, the whole file · a weak match · score 0.65 · post hoc, Riemannian distance, control animals, Holm, Cohen, G7
  5. [5] § Results › In silico modeling reveals complementary relationship between MePD GABA and UCN3 populations ↔ MePD_model/simulations/KNDyXMePD_stress.m, lines 3–56 · score 0.65 · neuronal activity, KNDy model, model parameter, couple, secretion, MePD
  6. [6] § Methods › Mini-endoscope-based calcium imaging › Restraint stress ↔ data_analysis/analysis_pipeline.ipynb, lines 2370–2403 · score 0.62 · GCaMP, calcium imaging, restraint stress, opsin, UCN3, GABA
  7. [7] § Methods › Modeling ↔ MePD_model/simulations/KNDyXMePDU.m, lines 2–61 · score 0.57 · model parameter, GABA efferent, couple, MePD, glutamate, KNDy
  8. [8] § Methods › Clustering analysis ↔ data_analysis/analysis_pipeline.ipynb, lines 1695–1714 · score 0.56 · linkage criterion, hierarchical clustering, agglomerative, activity
  9. [9] § Results › Selective optogenetic stimulation of UCN3 non-GABA neurons inhibits pulsatile LH secretion in female mice ↔ data_analysis/data/LH_profiling/figures/fig5J.m, the whole file · a weak match · score 0.54 · ConFon, LH IPI, bars, SEM, position, Figure 5
  10. [10] § Methods › Data analysis ↔ data_analysis/analysis_pipeline.ipynb, lines 495–514 · score 0.54 · power spectral density, GABA neurons, pre, max, signal, clustering
  11. [11] § Results › The MePD GABA neurons relay information from UCN3 neurons and inhibit the GnRH pulse generator ↔ data_analysis/data/LH_profiling/figures/fig5D.m, the whole file · a weak match · score 0.53 · ConFoff, LH IPI, vertical, bars, SEM, position
  12. [12] § Results › The MePD GABA neurons relay information from UCN3 neurons and inhibit the GnRH pulse generator ↔ data_analysis/data/LH_profiling/figures/fig5J.m, the whole file · a weak match · score 0.52 · ConFon, LH IPI, bars, SEM, position
  13. [13] § Methods › Mini-endoscope-based calcium imaging › Optogenetic stimulation of MePD UCN3 neurons ↔ data_analysis/analysis_pipeline.ipynb, lines 2370–2403 · score 0.52 · GCaMP, Calcium imaging, opsin, UCN3, GABA, stimulation
  14. [14] § Methods › Modeling ↔ MePD_model/simulations/combinedMePDv1.m, lines 43–73 · score 0.52 · GABA efferent neurons, GABA interneurons, opted, MePD, UCN3, model
  15. [15] § Results › Dynamic reconfiguration of MePD GABA neuronal networks during selective optogenetic stimulation of urocortin-3 (UCN3) neurons and restraint stress ↔ data_analysis/data/calcium_imaging/figures/permutation_test_AUC_data_fig1.ipynb, lines 9–33 · score 0.51 · stim stress, permutation, AUC, Cliff, bootstrap, delta
  16. [16] § Methods › Comparing pre-intervention and intervention cell-cell interaction strength ↔ data_analysis/calculate_riemannian_distances.ipynb, lines 150–172 · score 0.51 · Riemannian distance, correlation matrices, pre
  17. [17] § Results › Dynamic reconfiguration of MePD GABA neuronal networks during selective optogenetic stimulation of urocortin-3 (UCN3) neurons and restraint stress ↔ data_analysis/calculate_riemannian_distances.ipynb, lines 150–172 · score 0.51 · normalized Riemannian distance, correlation matrices, calcium

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 authors' code

Jupyter notebook · 2,816 lines · 104 KB · CC0-1.0 · 5 matches

  1. # %% [markdown]
  2. # # [Amygdala GABA Neurons: Gatekeepers of Stress and Reproduction in Female Mice ](https://www.biorxiv.org/content/10.1101/2025.01.06.631361v4). Analysis pipeline + Figures
  4. # %% [markdown]
  5. # # 🧰 Preparation Block for Analysis
  7. # %%
  8. # Check the current Python version for a stable connection
  9. !python --version
  11. # %% [markdown]
  12. # ### 📦 Install Required Packages
  14. # %%
  15. # Install the Shapely library for geometric operations
  16. !pip install shapely
  17. !pip install colabcode
  18. !pip install colorcet
  20. # Install scikit-posthocs for performing statistical post-hoc tests
  21. !pip install scikit_posthocs
  22. # Install statsmodels module for statistical modeling and hypothesis testing
  23. !python -m pip install statsmodels
  25. # %% [markdown]
  26. # ### 📥 Importing Required Packages
  28. # %%
  29. # OS interface for interacting with the file system
  30. import os
  32. # Core data manipulation and visualisation libraries
  33. import numpy as np # Numerical computing (arrays, math functions)
  34. import pandas as pd # Data structures & analysis (DataFrame, CSV I/O)
  35. import matplotlib.pyplot as plt # Plotting library for basic graphs
  36. import colorcet as cc # Colour maps for visualisation
  37. import seaborn as sns # Statistical data visualisation, builds on Matplotlib
  38. from pathlib import Path # Object-oriented filesystem paths
  40. # Advanced plotting tools from Matplotlib
  41. from matplotlib import gridspec # For custom subplot layouts
  42. from matplotlib.cm import ScalarMappable # For colour mapping in plots
  44. # Signal processing & statistical functions from SciPy
  45. from scipy.signal import correlate, periodogram # Signal similarity & frequency analysis
  46. from scipy.cluster.hierarchy import linkage, dendrogram # Hierarchical clustering & visual trees
  47. from scipy.spatial.distance import pdist, squareform # Pairwise distance computations
  48. from scipy.stats import mannwhitneyu, kstest, kruskal # Non-parametric tests
  49. from scipy.stats import wilcoxon # Wilcoxon signed-rank test
  50. from scipy.linalg import logm, sqrtm # Matrix operations for Riemannian geometry
  51. from scipy.stats import t
  53. # Machine learning tools from Scikit-learn
  54. from sklearn.preprocessing import scale, minmax_scale # Feature scaling methods
  55. from sklearn.cluster import KMeans, AgglomerativeClustering # Clustering algorithms
  56. from sklearn.metrics import silhouette_score # Clustering quality metric
  58. import statsmodels.formula.api as smf # Statistical modeling using formulas
  60. # Utilities
  61. from itertools import product # Cartesian product generator (combinatorics)
  63. # Shapes for visualisation
  64. from matplotlib.patches import Ellipse, Polygon # Custom shapes for plots
  65. from shapely.geometry import Polygon as ShapelyPolygon # Geometric operations
  67. # Post-hoc statistical tests
  68. from scikit_posthocs import posthoc_dunn # Dunn’s test for multiple comparisons after Kruskal-Wallis
  70. # Riemannian distance function
  71. from pyriemann.utils.distance import distance_riemann
  73. # ✅ Confirmation
  74. print("Imports done ✅!")
  76. # %% [markdown]
  77. # These two sections below are optional. They are only required for consistency of the style of the figures.
  79. # %%
  80. # Set the default appearance for Seaborn plots to enhance readability and aesthetics
  81. sns.set_theme(
  82. context='notebook', # Context presets size and scaling for notebook display
  83. style='whitegrid', # Background grid style for better visual separation
  84. palette='deep', # colour palette for plots (rich, vivid colours)
  85. font='DejaVu Sans', # Font for labels and titles
  86. font_scale=1.5, # Increase font size for better legibility
  87. rc={ # Override specific matplotlib parameters
  88. 'lines.linewidth': 2, # Thicker lines for better visibility
  89. 'lines.markersize': 8, # Larger markers for clarity
  90. 'xtick.labelsize': 15, # Font size for x-axis tick labels
  91. 'ytick.labelsize': 15 # Font size for y-axis tick labels
  92. }
  93. )
  95. # %%
  96. # Customise Matplotlib settings for consistent plot styling
  97. rc_params = {
  98. 'axes.grid': False, # Turn off background gridlines
  99. 'font.family': 'DejaVu Sans', # Set global font for plot text
  100. 'font.size': 15, # Base font size for all elements
  101. 'lines.linewidth': 2, # Make plot lines thicker and easier to read
  102. 'lines.markersize': 8, # Larger markers for better visibility
  103. 'xtick.major.size': 3, # Length of major ticks on x-axis
  104. 'ytick.major.size': 3, # Length of major ticks on y-axis
  105. 'xtick.labelsize': 15, # Font size for x-axis tick labels
  106. 'ytick.labelsize': 15, # Font size for y-axis tick labels
  107. 'xtick.bottom': True, # Show ticks on bottom of x-axis
  108. 'xtick.top': False, # Hide ticks on top of x-axis
  109. 'ytick.left': True, # Show ticks on left of y-axis
  110. 'ytick.right': False, # Hide ticks on right of y-axis
  111. 'axes.edgecolor': 'grey' # Use grey colour for plot borders
  112. }
  114. # Apply each setting to Matplotlib's runtime configuration
  115. for param, value in rc_params.items():
  116. plt.rcParams[param] = value
  118. # %% [markdown]
  119. # ### 🛠 Utility Functions for Your Analysis
  121. # %% [markdown]
  122. # Several functions have been defined here to avoid repeating the code; however, in parts, I have repeated these steps with appropriate comments for better readability.
  124. # %%
  125. def load_data(stim_csv_fname):
  126. """
  127. Load stimulation data from a CSV file and extract relevant cell trace information.
  129. Parameters:
  130. stim_csv_fname (str): Path to the CSV file containing stimulation data.
  132. Returns:
  133. tuple:
  134. - data_stim (pd.DataFrame): Entire transposed dataset including metadata.
  135. - cells_traces_tot_old (pd.DataFrame): Data subset excluding metadata (typically fluorescence traces).
  136. """
  138. # Load data from CSV and transpose it (rows become columns and vice versa)
  139. data_stim = pd.read_csv(stim_csv_fname).T
  141. # Strip leading/trailing whitespace from index labels to standardise them
  142. data_stim.index = data_stim.index.str.strip()
  144. # Extract cell trace data by skipping the first row and column (assumed to be metadata)
  145. cells_traces_tot_old = data_stim.iloc[1:, 1:]
  147. # Return both the full dataset and the cleaned subset for downstream use
  148. return data_stim, cells_traces_tot_old
  150. # %%
  151. def process_data(original_data, processed_data, stim_onset=1204, stim_offset=2404, stimulus_freq=0.1, kmeans_clusters=2, random_state=11111):
  152. """
  153. Process fluorescence stimulation data to identify and remove noisy/resonating cells
  154. using frequency-based filtering and K-means clustering.
  156. Parameters:
  157. original_data (pd.DataFrame): Raw stimulation data including metadata.
  158. processed_data (pd.DataFrame): Data with metadata removed (fluorescence traces).
  159. stim_onset (int): Index of stimulation onset.
  160. stim_offset (int): Index of stimulation offset.
  161. stimulus_freq (float): Expected stimulation frequency in Hz (default is 0.1 Hz).
  162. kmeans_clusters (int): Number of clusters used in K-means (default is 2).
  163. random_state (int): Seed for reproducible clustering (default is 11111).
  165. Returns:
  166. tuple:
  167. - data_stim_tot (np.ndarray): Cleaned fluorescence data with noisy cells removed.
  168. - n_cells (int): Number of remaining, reliable cells.
  169. """
  171. # Get number of cells (rows) from processed data
  172. total_cells_stim, _ = processed_data.shape
  174. # Extract time data from original_data and calculate sampling frequency
  175. ts_stim_tot = np.asarray(original_data.iloc[0, 1:], dtype=float)
  176. sec_per_frame = ts_stim_tot[1]
  177. fs = 1 / sec_per_frame # Sampling frequency in Hz
  179. # Convert processed data (traces) to NumPy array for efficient computation
  180. data_stim_artifact = np.asarray(processed_data, dtype=float)
  182. # Identify resonating cells via frequency analysis (Periodogram)
  183. labels_removed_PSD = []
  184. for i in range(total_cells_stim):
  185. freq, Pxx_den = periodogram(data_stim_artifact[i, stim_onset:stim_offset], fs)
  186. Pxx_den = Pxx_den / max(Pxx_den) # Normalise power
  187. freq_bin = round(stimulus_freq * (stim_offset - stim_onset) / fs)
  188. if Pxx_den[freq_bin] > 0.37:
  189. labels_removed_PSD.append(i) # Mark as resonating
  191. # Create binary mask of valid vs. resonating cells
  192. labels_artifact = np.array([0 if i in labels_removed_PSD else 1 for i in range(total_cells_stim)])
  194. # Apply K-means clustering to activity during stimulus period
  195. km = KMeans(n_clusters=kmeans_clusters, random_state=random_state, n_init=10, max_iter=100)
  196. clustered = scale(data_stim_artifact[:, stim_onset:stim_offset], axis=0)
  197. km.fit(clustered)
  198. labels_artifact_KM = km.labels_
  200. # Identify which cluster is smaller (assumed artifact) and remove it
  201. cluster_removed = 0 if np.sum(labels_artifact_KM == 0) < np.sum(labels_artifact_KM == 1) else 1
  202. for i in range(len(data_stim_artifact)):
  203. if labels_artifact[i] == 1 and labels_artifact_KM[i] == cluster_removed:
  204. labels_artifact[i] = 0 # Remove cell
  206. # Extract indices of removed cells and drop them
  207. labels_removed = np.where(labels_artifact == 0)[0]
  208. removed_row_labels = processed_data.index[labels_removed]
  209. cells_traces_stim_tot = processed_data.drop(removed_row_labels, axis=0, inplace=False)
  211. # Convert final cleaned DataFrame to NumPy array
  212. data_stim_tot = np.asarray(cells_traces_stim_tot, dtype=float)
  214. # Return cleaned data and count of remaining cells
  215. n_cells = len(cells_traces_stim_tot)
  216. return data_stim_tot, n_cells
  218. # %%
  219. def min_max_norm(data):
  220. # Convert input to NumPy array to ensure consistent operations
  221. data = np.asarray(data)
  223. # Compute per-trace (row-wise) minimum values
  224. min_vals = np.min(data, axis=1, keepdims=True)
  226. # Compute per-trace (row-wise) maximum values
  227. max_vals = np.max(data, axis=1, keepdims=True)
  229. # Scale each trace independently to the range [0, 1]
  230. return (data - min_vals) / (max_vals - min_vals)
  233. # %%
  234. # Convert p-value to custom significance symbols for visualisation
  235. def pval_to_star(pvalue):
  236. """
  237. Converts a p-value into a string representing significance levels using hash marks and daggers.
  239. Parameters:
  240. pvalue (float): The p-value to evaluate.
  242. Returns:
  243. str: A string indicating statistical significance:
  244. '##' → Highly significant (p ≤ 0.001)
  245. '#' → Significant (p ≤ 0.01)
  246. '†' → Marginally significant (p ≤ 0.05)
  247. 'ns' → Not significant
  248. """
  249. if pvalue <= 0.0001:
  250. return '##' # Extremely significant
  251. elif pvalue <= 0.001:
  252. return '##' # Highly significant
  253. elif pvalue <= 0.01:
  254. return '#' # Significant
  255. elif pvalue <= 0.05:
  256. return '†' # Marginally significant
  257. else:
  258. return 'ns' # Not significant
  260. # %%
  261. # Compute a metric based on signed lagged cross-correlation (SLxCorr)
  262. def signed_lagged_corr(x, y, demean=True):
  263. """
  264. Computes the signed lagged cross-correlation between two signals and returns the
  265. correlation value at the lag with the maximum absolute correlation, along with that lag.
  267. Parameters:
  268. x (numpy.ndarray): First input signal (1D).
  269. y (numpy.ndarray): Second input signal (1D).
  270. demean (bool): If True, subtract the mean from each signal before correlation.
  272. Returns:
  273. tuple:
  274. - float: Signed, normalised cross-correlation at the peak absolute correlation.
  275. - int: Lag (in samples) corresponding to that peak. Positive lag means y lags x.
  276. """
  277. # Ensure 1D arrays
  278. x = np.asarray(x).ravel()
  279. y = np.asarray(y).ravel()
  281. # Optionally remove DC offset so correlation focuses on shape, not level
  282. if demean:
  283. x = x - x.mean()
  284. y = y - y.mean()
  286. # Guard against zero-variance inputs
  287. norm = np.linalg.norm(x) * np.linalg.norm(y)
  288. if norm == 0:
  289. raise ValueError("Normalisation factor is zero. Ensure signals have non-zero variance.")
  291. # Full cross-correlation; zero lag is at index len(y) - 1 for correlate(x, y)
  292. cross_corr = correlate(x, y, mode='full')
  293. scc = cross_corr / norm # Normalised cross-correlation
  294. abs_scc = np.abs(scc)
  296. # Peak absolute correlation and its lag
  297. max_idx = int(np.argmax(abs_scc))
  298. lag = max_idx - (len(y) - 1) # Convert index to lag (centered at zero lag)
  299. signed_correlation = float(scc[max_idx]) # Keep the sign at the peak
  301. return signed_correlation, lag
  303. # %%
  304. import warnings
  305. warnings.filterwarnings("ignore", category=RuntimeWarning)
  306. # run the risky operations here
  308. # Generate mean and standard deviation of Riemannian distances for SPD (correlation) matrices of sizes p=5–80
  309. def metropolis_hastings(p, component_index, num_samples=2000, burn_in=200, proposal_std=0.01):
  310. """
  311. Generate random correlation matrix components using Metropolis–Hastings sampling on the unit sphere.
  313. Parameters:
  314. p (int): Dimension of the full correlation matrix.
  315. component_index (int): Index of the component being generated (1-based).
  316. num_samples (int): Number of samples to retain after burn-in.
  317. burn_in (int): Number of initial samples to discard to allow convergence.
  318. proposal_std (float): Standard deviation of the Gaussian proposal noise.
  320. Returns:
  321. numpy.ndarray: Array of shape (num_samples, dim) containing generated unit vectors.
  322. """
  323. dim = p - component_index + 1 # Effective dimension of the current component
  324. samples = np.zeros((burn_in + num_samples + 1, dim))
  326. # Initialise with a random unit vector; ensure first component is non-negative
  327. current = np.random.randn(1, dim)
  328. current[0, 0] = np.abs(current[0, 0])
  329. current = current / np.linalg.norm(current)
  331. # Run Metropolis–Hastings chain
  332. for t in range(burn_in + num_samples + 1):
  333. # Propose a new vector with small Gaussian perturbation
  334. proposal = current + np.random.normal(scale=proposal_std, size=(1, dim))
  335. proposal = proposal / np.linalg.norm(proposal) # Renormalise to unit length
  337. # Compute acceptance ratio (ensuring positivity of first component)
  338. delta = np.random.uniform(0, 1)
  339. acceptance_ratio = (proposal[0, 0] / current[0, 0])**component_index if proposal[0, 0] >= 0 else 0
  341. # Accept or reject proposal
  342. if delta < acceptance_ratio:
  343. current = proposal
  345. samples[t, :] = current.flatten()
  347. # Return only the post-burn-in samples
  348. return samples[burn_in + 1:, :]
  351. def random_correlation_matrix(p):
  352. """
  353. Construct a random correlation matrix of size p × p using Metropolis–Hastings sampling.
  355. Each row of the upper-triangular matrix U is sampled sequentially, and the
  356. final correlation matrix is obtained as C = U Uᵀ, normalised to have unit diagonal.
  358. Parameters:
  359. p (int): Dimension of the correlation matrix.
  361. Returns:
  362. numpy.ndarray: Random correlation matrix of size p × p.
  363. """
  364. U = np.zeros((p, p))
  365. for i in range(p):
  366. # For each row i, sample the corresponding unit vector and use the last draw
  367. U[i, i:] = metropolis_hastings(p, i + 1)[-1, :]
  369. # Compute symmetric product and normalise diagonal entries to 1
  370. C = U @ U.T
  371. D = np.sqrt(np.diag(C))
  372. correlation_matrix = C / np.outer(D, D)
  373. return correlation_matrix
  376. def distance_riemann(A, B):
  377. """
  378. Compute the affine-invariant Riemannian distance between two SPD matrices A and B.
  380. Formula:
  381. d(A, B) = || log( A^{-1/2} B A^{-1/2} ) ||_F
  383. Parameters:
  384. A (numpy.ndarray): First SPD matrix.
  385. B (numpy.ndarray): Second SPD matrix.
  387. Returns:
  388. float: The affine-invariant Riemannian distance between the two matrices.
  389. """
  390. sqrt_A = sqrtm(A)
  391. inv_sqrt_A = np.linalg.inv(sqrt_A)
  392. C = inv_sqrt_A @ B @ inv_sqrt_A
  393. log_C = logm(C)
  394. return np.linalg.norm(log_C, 'fro')
  397. def expected_distance(matrix_size, num_pairs=100):
  398. """
  399. Estimate the expected Riemannian distance (mean ± std) between pairs of
  400. random correlation matrices of a given size.
  402. Parameters:
  403. matrix_size (int): Dimension of the correlation matrices.
  404. num_pairs (int): Number of pairs of matrices to generate.
  406. Returns:
  407. tuple:
  408. - float: Mean of the estimated distances.
  409. - float: Standard deviation of the estimated distances.
  410. """
  411. distances = []
  412. for _ in range(num_pairs):
  413. A = random_correlation_matrix(matrix_size)
  414. B = random_correlation_matrix(matrix_size)
  415. d = distance_riemann(A, B)
  416. distances.append(d)
  418. return np.mean(distances), np.std(distances)
  422. # %%
  423. def compute_AUC(data_tot, stim_onset, stim_offset):
  424. """
  425. Compute the change in area under the curve (AUC) of ΔF/F signals before and during stimulation.
  427. Parameters:
  428. data_tot (numpy.ndarray): 2D array of calcium fluorescence signals
  429. (shape: n_cells × n_timepoints).
  430. stim_onset (int): Frame index marking the start of stimulation.
  431. stim_offset (int): Frame index marking the end of stimulation.
  433. Notes:
  434. - ΔF/F (delta F over F) is computed relative to the mean baseline fluorescence
  435. during the pre-stimulation period.
  436. - AUC is computed for several time windows (30, 60, 90, 120 seconds)
  437. within the stimulation period.
  438. - Sampling rate is assumed to be 10 Hz (hence `time_point * 10`).
  439. """
  441. AUC_time_windows = [30, 60, 90, stim_offset/10] # time windows (in seconds) for which AUC is computed
  443. # --- Pre-stimulation ΔF/F computation ---
  444. # Compute mean baseline for each cell during pre-stimulation period
  445. baseline_pre = np.mean(data_tot[:, :stim_onset], axis=1)
  446. # Compute ΔF/F trace for pre-stimulation
  447. dff_pre = (data_tot[:, :stim_onset] - baseline_pre[:, np.newaxis]) / baseline_pre[:, np.newaxis]
  449. # --- Stimulation ΔF/F computation ---
  450. # Compute mean baseline again (same as above, to ensure comparable scaling)
  451. baseline_stim = np.mean(data_tot[:, :stim_onset], axis=1)
  452. # Compute ΔF/F trace during stimulation
  453. dff_stim = (data_tot[:, stim_onset:stim_offset] - baseline_stim[:, np.newaxis]) / baseline_stim[:, np.newaxis]
  455. # --- Compute AUC differences for each time window ---
  456. for AUC_window in AUC_time_windows:
  457. # Compute average AUC for pre-stimulation (normalised by duration)
  458. area_1 = np.trapezoid(dff_pre) / (stim_onset / 10)
  460. # Compute average AUC during stimulation (up to specified duration)
  461. # Sampling rate assumed 10 Hz → convert seconds to frame index.
  462. area_2 = np.trapezoid(dff_stim[:, :AUC_window * 10]) / AUC_window
  464. # Output median difference in AUC between stimulation and pre-stimulation
  465. print(f'AUC for {AUC_window} s is {np.median(area_2) - np.median(area_1)}.')
  467. # %% [markdown]
  468. # ### 📁 Setting File Paths for Figures and Data
  470. # %%
  471. # --- Fig. 1: Stimulation and Stress Data ---
  472. stim_fig1_csv_fname = './data/calcium_imaging/UCN3_stimulation/G7 20240201 10hz 223sti.csv'
  473. stress_fig1_csv_fname = './data/calcium_imaging/restrain_stress/G28 333restraint 20240630.csv'
  475. # --- Fig. 2: Stimulation and Stress Data ---
  476. stim_fig2_csv_fname = './data/calcium_imaging/UCN3_stimulation/G28 20240709 10hz 223sti.csv'
  477. stress_fig2_csv_fname = './data/calcium_imaging/restrain_stress/G9 20231218 334restraint.csv'
  479. # --- Figs. 1 & 2: P-values and Cell Percentages ---
  480. hier_clustering_cell_perc_fname = './data/cell_percentage_hier_corr 10Hz.csv'
  483. if not Path(stim_fig1_csv_fname).is_file():
  484. print(f"File {stim_fig1_csv_fname} does not exist")
  485. if not Path(stress_fig1_csv_fname).is_file():
  486. print(f"File {stim_fig1_csv_fname} does not exist")
  487. if not Path(stim_fig2_csv_fname).is_file():
  488. print(f"File {stim_fig1_csv_fname} does not exist")
  489. if not Path(stress_fig2_csv_fname).is_file():
  490. print(f"File {stim_fig1_csv_fname} does not exist")
  491. if not Path(hier_clustering_cell_perc_fname).is_file():
  492. print(f"File {hier_clustering_cell_perc_fname} does not exist")
  495. # %%
  496. # --- Time Windows (in index units) ---
  497. stim_onset = 1204
  498. stim_offset = 2404
  499. stress_onset = 1804
  500. stress_offset = 3604
  502. # --- Signal Analysis ---
  503. PSD_thr = 0.37 # Power Spectral Density threshold
  505. # --- Clustering & Bootstrapping ---
  506. max_no_clusters = 4
  507. n_resampling = 100
  508. Ks_1 = 2 # Pre-stimulation clusters
  509. Ks_2 = 2 # Stimulation clusters
  511. # --- GABA Neuron colours ---
  512. eff_green = (63/255, 128/255, 0) # GABA efferent
  513. inter_green = (0, 1, 0) # GABA interneuron
  514. GABA_greens = [eff_green, inter_green]
  516. # %% [markdown]
  517. # # 🧪 Full Workflow: G7 Stimulation (Figures 1 & S1)
  519. # %% [markdown]
  520. # This is an example stimulation trial for figure 1 from animal G7.
  522. # %%
  523. # Load and transpose the data
  524. data_fig1_stim = pd.read_csv(stim_fig1_csv_fname).T
  526. # Clean up row labels
  527. data_fig1_stim.index = data_fig1_stim.index.str.strip()
  529. # Preview the last 5 cells
  530. data_fig1_stim.tail(5)
  532. # %% [markdown]
  533. # ## ✂️ Data Slicing for G7 Stimulation
  535. # %%
  536. # Extract signal data (excluding metadata)
  537. cells_traces_tot_old = data_fig1_stim.iloc[1:, 1:]
  539. # Get shape info
  540. total_cells_fig1_stim, data_fig1_stim_length_tot = cells_traces_tot_old.shape
  542. # Preview bottom rows
  543. cells_traces_tot_old.tail()
  545. # %%
  546. # Extract time axis from first row (excluding first column)
  547. ts_fig1_stim_tot = np.asarray(data_fig1_stim.iloc[0, 1:], dtype=float)
  549. # Time between frames (first timestamp is 0)
  550. sec_per_frame = ts_fig1_stim_tot[1]
  552. # Sampling frequency (Hz)
  553. fs = 1 / sec_per_frame
  555. # %%
  556. # Convert fluorescence traces to NumPy array
  557. data_fig1_stim_artifact = np.asarray(cells_traces_tot_old, dtype=float)
  559. # %% [markdown]
  560. # ## 🔍 Finding UCN3 Cells
  562. # %% [markdown]
  563. # ### 🌀 1. Power Spectral Density (PSD) Filtering
  565. # %%
  566. # Define the stimulus frequency and initialise a list to store labels of removed cells based on PSD threshold
  567. stimulus_freq = 0.1 # Stimulus frequency in Hz
  568. labels_removed_PSD = [] # List to store indices of cells flagged as responsive
  570. # Loop through each cell to calculate the power spectral density (PSD)
  571. for i in range(total_cells_fig1_stim):
  572. # Compute the periodogram (PSD) for the cell's signal during the stimulation period
  573. freq, Pxx_den = periodogram(
  574. data_fig1_stim_artifact[i, stim_onset:stim_offset], fs)
  576. # Normalise the PSD by its maximum value to scale between 0 and 1
  577. Pxx_den = Pxx_den / max(Pxx_den)
  579. # Calculate the index corresponding to the stimulus frequency
  580. # This assumes uniform frequency spacing in the PSD output
  581. stimulus_freq_idx = round(stimulus_freq * (stim_offset - stim_onset) / fs)
  583. # Check if the normalised PSD at the stimulus frequency exceeds the threshold
  584. if Pxx_den[stimulus_freq_idx] > PSD_thr:
  585. # If so, mark the cell as responsive by adding its index to the list
  586. labels_removed_PSD.append(i)
  588. # Print the indices of cells identified as responsive based on PSD threshold
  589. print("Indices of cells removed based on PSD threshold:", labels_removed_PSD)
  591. # %%
  592. # Create binary mask to flag artifact-affected cells
  593. # - 0: cell identified as artifact (high PSD at stimulus frequency)
  594. # - 1: cell retained for further analysis
  595. labels_artifact = np.asarray(
  596. [0 if _ in labels_removed_PSD else 1 for _ in range(total_cells_fig1_stim)]
  597. )
  599. # Optional: Print the binary mask to verify which cells were flagged
  600. # 0 = artifact (failed PSD threshold), 1 = valid (passed PSD threshold)
  601. print(labels_artifact)
  603. # %% [markdown]
  604. # ### 🔢 2. K-Means Clustering
  606. # %%
  607. # Apply scikit-learn K-means clustering to sort activity traces
  609. # Step 1: Normalise the activity data across time (columns) to ensure fair clustering
  610. scaled_data = scale(data_fig1_stim_artifact[:, stim_onset:stim_offset], axis=0)
  612. # Step 2: Set the number of clusters (2 groups: responsive vs non-responsive)
  613. Ks = 2
  615. # Step 3: Initialise and run K-means clustering
  616. # - random_state ensures reproducibility
  617. # - n_init = 10 means the algorithm runs 10 times with different centroid seeds
  618. # - max_iter = 100 limits the number of iterations per run
  619. kmeans = KMeans(n_clusters=Ks, random_state=11111, n_init=10, max_iter=100)
  621. # Step 4: Fit the model to the scaled data and assign cluster labels to each cell
  622. labels_artifact_KM = kmeans.fit_predict(scaled_data)
  624. # Step 5: Sort cell indices based on their assigned cluster labels
  625. # This helps organise or visualise cells by their activity patterns
  626. i_labels_artifact = np.argsort(labels_artifact_KM)
  628. # Step 6: Display the sorted indices for inspection or downstream analysis
  629. print(f"Sorted indices of cells by cluster ID:\n{i_labels_artifact}")
  631. # %%
  632. # Determine which cluster to remove based on relative size
  634. # Step 1: Count the number of cells in each cluster
  635. size_cluster_0 = len(data_fig1_stim_artifact[labels_artifact_KM == 0])
  636. size_cluster_1 = len(data_fig1_stim_artifact[labels_artifact_KM == 1])
  638. # Step 2: Compare cluster sizes to identify the smaller one
  639. # The assumption here is that the smaller cluster likely represents artifacts or noise
  640. if size_cluster_0 < size_cluster_1:
  641. cluster_removed = 0 # Smaller cluster is flagged for removal
  642. cluster_kept = 1 # Larger cluster is retained
  643. print('Cluster #0 (smaller) is removed.')
  644. elif size_cluster_1 != 0:
  645. cluster_removed = 1
  646. cluster_kept = 0
  647. print('Cluster #1 (smaller) is removed.')
  648. else:
  649. # Edge case: one cluster has zero members, so no valid removal decision
  650. print("No valid cluster for removal.")
  652. # %%
  653. # Update artifact labels based on K-means clustering results
  655. # Step 1: Check if the cluster marked for removal contains any cells
  656. if len(data_fig1_stim_artifact[labels_artifact_KM == cluster_removed]) != 0:
  657. # Step 2: Iterate through all cells
  658. for i in range(len(data_fig1_stim_artifact)):
  659. # Only update cells that were previously marked as valid (1)
  660. # and now belong to the cluster flagged for removal
  661. if labels_artifact[i] == 1 and labels_artifact_KM[i] == cluster_removed:
  662. labels_artifact[i] = 0 # Mark as artifact
  664. # Step 3: Explicitly define which cluster was removed and which was kept
  665. cluster_removed = 0
  666. cluster_kept = 1
  668. # Step 4: Extract indices of cells now marked as artifacts
  669. labels_removed = np.where(labels_artifact == cluster_removed)[0]
  671. # Optional: print summary for verification
  672. print(f"Number of (UCN) cells removed: {len(labels_removed)}")
  673. print(f"Indices of removed cells: {labels_removed}")
  675. # %%
  676. # Remove resonated (artifact) cells from fluorescence intensity traces
  678. # Step 1: Get the row labels (e.g., cell IDs) corresponding to removed indices
  679. removed_row_labels = cells_traces_tot_old.index[labels_removed]
  681. # Step 2: Drop the identified artifact cells from the original DataFrame
  682. # - axis=0 means rows are being removed
  683. # - inplace=False ensures the original DataFrame remains unchanged
  684. cells_traces_fig1_stim_tot = cells_traces_tot_old.drop(
  685. removed_row_labels, axis=0, inplace=False)
  687. # Step 3: Count the number of remaining (valid) cells after removal
  688. n_cells_fig1_stim = len(cells_traces_fig1_stim_tot)
  690. # Optional: print summary
  691. print(f"Number of GABA cells: {n_cells_fig1_stim}")
  693. # %% [markdown]
  694. # ### ✂️ Data Slicing
  696. # %%
  697. # Convert filtered fluorescence traces from DataFrame to NumPy array
  698. data_fig1_stim_tot = cells_traces_fig1_stim_tot.to_numpy(dtype=float)
  700. # Slice the data into temporal segments for analysis
  701. # Pre-stimulation period: from start up to stimulus onset
  702. data_fig1_stim_1 = data_fig1_stim_tot[:, :stim_onset]
  704. # Stimulation period: from stimulus onset to stimulus offset
  705. data_fig1_stim_2 = data_fig1_stim_tot[:, stim_onset:stim_offset]
  707. # Optional: Print shapes to verify correct slicing
  708. print(f"Total data shape: {data_fig1_stim_tot.shape}")
  709. print(f"Pre-stimulation data shape: {data_fig1_stim_1.shape}")
  710. print(f"Stimulation data shape: {data_fig1_stim_2.shape}")
  712. # %%
  713. # Rescale pre-stimulation traces using min-max normalisation (row-wise)
  714. normalised_data_fig1_stim_1 = minmax_scale(data_fig1_stim_1, axis=1)
  716. # Optional: Inspect the first few normalised traces to verify the transformation
  717. print("Normalised Pre-stimulation Data (first few rows):")
  718. print(normalised_data_fig1_stim_1[:5])
  720. # %%
  721. # Rescale stimulation-period traces using min-max normalisation (row-wise)
  722. normalised_data_fig1_stim_2 = minmax_scale(data_fig1_stim_2, axis=1)
  724. # Optional: Inspect the first few normalised stimulation traces
  725. print("Normalised Stimulation Data (first few rows):")
  726. print(normalised_data_fig1_stim_2[:5])
  728. # %% [markdown]
  729. # ### 📈 corrcoef: Standard Pearson Correlation Coefficient
  731. # %%
  732. # Compute the correlation matrix for pre-stimulation traces
  733. max_corr_fig1_stim_1 = np.corrcoef(normalised_data_fig1_stim_1, rowvar=True)
  735. # Compute the correlation matrix for stimulation-period traces
  736. max_corr_fig1_stim_2 = np.corrcoef(normalised_data_fig1_stim_2, rowvar=True)
  738. # %% [markdown]
  739. # ## 🖼️ Figure S2: Visualisation of Trace Filtering and Clustering
  741. # %%
  742. # Set up the figure canvas with 4 vertical subplots
  743. plt.figure(figsize=(15, 14))
  745. # ─────────────────────────────────────────────────────────────
  746. # 🔹 Plot 1: Original Traces (before any filtering)
  747. ax1 = plt.subplot(4, 1, 1)
  749. # Normalise raw traces for visualisation
  750. scaled_data = minmax_scale(data_fig1_stim_artifact[:, :stim_offset], axis=1)
  752. # Display as heatmap
  753. plt.imshow(scaled_data, cmap="viridis", interpolation='none')
  755. # Mark stimulus onset with a vertical white line
  756. ax1.plot([1204, 1204], [0, total_cells_fig1_stim-1], '-w')
  758. # Configure plot aesthetics
  759. ax1.set_aspect('auto')
  760. ax1.set_xlabel('Time [sec]')
  761. ax1.set_ylabel('Cell ID')
  762. ax1.invert_yaxis()
  763. ax1.set_title('Original Traces')
  764. ax1.set_xticks([0, 500, 1000, 1500, 2000])
  765. ax1.set_xticklabels([0, 50, 100, 150, 200])
  767. # Add colourbar for ΔF/F₀ intensity
  768. cmap = plt.get_cmap("viridis")
  769. sm = ScalarMappable(cmap=cmap)
  770. sm.set_clim(vmin=0, vmax=1)
  771. cb = plt.colorbar(cmap=cmap, fraction=0.04, pad=0.03)
  772. cb.set_label('$\Delta$ F$/$F$_0$ (scaled)')
  773. cb.set_ticks([0, 0.5, 1])
  775. # ─────────────────────────────────────────────────────────────
  776. # 🔹 Plot 2: Power Spectral Density (PSD) at 0.1 Hz
  777. ax2 = plt.subplot(4, 1, 2)
  779. # Loop through each cell to compute and plot PSD
  780. for i in range(total_cells_fig1_stim):
  781. freq, Pxx_den = periodogram(data_fig1_stim_artifact[i, stim_onset:stim_offset], fs)
  782. Pxx_den = Pxx_den / max(Pxx_den) # Normalise PSD
  783. stimulus_freq_PSD = Pxx_den[round(stimulus_freq * (stim_offset - stim_onset) / fs)]
  785. # Mark cells based on PSD threshold
  786. if stimulus_freq_PSD < PSD_thr:
  787. labels_removed_PSD.append(i)
  788. ax2.scatter(i, stimulus_freq_PSD, marker='.', s=300, c='black') # Below threshold
  789. else:
  790. ax2.scatter(i, stimulus_freq_PSD, marker='.', s=300, c='red') # Above threshold
  792. # Add horizontal threshold line
  793. ax2.plot([-1, total_cells_fig1_stim], [0.37, 0.37], '--', c='gray')
  794. ax2.set_xlim([-1, total_cells_fig1_stim])
  795. ax2.set_ylim([-0.05, 1.1])
  796. ax2.set_xlabel('Cell ID')
  797. ax2.set_ylabel('PSD')
  798. ax2.set_title('Power Spectral Density (PSD) at 0.1 Hz')
  799. ax2.annotate(0.37, (10, 0.4), ha='center')
  801. # ─────────────────────────────────────────────────────────────
  802. # 🔹 Plot 3: Traces Sorted by K-means Clustering
  803. ax3 = plt.subplot(4, 1, 3)
  805. # Display sorted traces based on clustering labels
  806. plt.imshow(minmax_scale(data_fig1_stim_artifact[i_labels_artifact, :stim_offset], axis=1), cmap="viridis", interpolation='none')
  808. # Mark stimulus onset
  809. ax3.plot([1204, 1204], [0, total_cells_fig1_stim-1], '-w')
  811. # Configure plot aesthetics
  812. ax3.set_aspect('auto')
  813. ax3.set_xlabel('Time [sec]')
  814. ax3.set_ylabel('Cell ID')
  815. ax3.invert_yaxis()
  816. ax3.set_title('K-means Sorted Traces')
  817. ax3.set_xticks([0, 500, 1000, 1500, 2000])
  818. ax3.set_xticklabels([0, 50, 100, 150, 200])
  820. # Add colourbar
  821. cmap = plt.get_cmap("viridis")
  822. sm = ScalarMappable(cmap=cmap)
  823. sm.set_clim(vmin=0, vmax=1)
  824. cb = plt.colorbar(cmap=cmap, fraction=0.04, pad=0.03)
  825. cb.set_label('$\Delta$ F$/$F$_0$ (scaled)')
  826. cb.set_ticks([0, 0.5, 1])
  828. # ─────────────────────────────────────────────────────────────
  829. # 🔹 Plot 4: Purely GABA Traces (after filtering)
  830. ax4 = plt.subplot(4, 1, 4)
  832. # Display final cleaned traces
  833. plt.imshow(minmax_scale(data_fig1_stim_tot[:, :stim_offset], axis=1), cmap="viridis", interpolation='none')
  835. # Mark stimulus onset
  836. ax4.plot([1204, 1204], [0, n_cells_fig1_stim-1], '-w')
  838. # Configure plot aesthetics
  839. ax4.set_aspect('auto')
  840. ax4.set_xlabel('Time [sec]')
  841. ax4.set_ylabel('Cell ID')
  842. ax4.invert_yaxis()
  843. ax4.set_title('Purely GABA Traces')
  844. ax4.set_xticks([0, 500, 1000, 1500, 2000])
  845. ax4.set_xticklabels([0, 50, 100, 150, 200])
  847. # Add colourbar
  848. cmap = plt.get_cmap("viridis")
  849. sm = ScalarMappable(cmap=cmap)
  850. sm.set_clim(vmin=0, vmax=1)
  851. cb = plt.colorbar(cmap=cmap, fraction=0.04, pad=0.03)
  852. cb.set_label('$\Delta$ F$/$F$_0$ (scaled)')
  853. cb.set_ticks([0, 0.5, 1])
  855. # ─────────────────────────────────────────────────────────────
  856. # Final layout adjustments and export
  857. plt.tight_layout()
  859. # Optional: Save figure to disk
  860. # plt.savefig("FigS1_revised.svg", dpi=600)
  862. # Display the figure
  863. plt.show()
  865. # %%
  866. # Visualise a few example stimulus traces after normalisation
  867. X = minmax_scale(data_fig1_stim_artifact[:, :stim_offset], axis=1)
  869. plt.figure(figsize=(10, 6))
  870. for i in range(3, 5): # plot a small subset of neurons
  871. plt.plot(X[i, :], linewidth=1)
  873. plt.xlim(0, stim_offset)
  874. plt.xlabel("Time")
  875. plt.ylabel("Normalised activity")
  876. plt.title("Stimulus traces (neurons 3–7)")
  877. plt.show()
  880. # %% [markdown]
  881. # # 🧪 Full Workflow: G28 Stress (Figure 1)
  883. # %% [markdown]
  884. # This is a example restraint stress trial for figure 1 from animal G28.
  886. # %%
  887. # Load and transpose the data
  888. data_fig1_stress = pd.read_csv(stress_fig1_csv_fname).T
  890. # Clean up row labels
  891. data_fig1_stress.index = data_fig1_stress.index.str.strip()
  893. # Preview the last 6 cells
  894. data_fig1_stress.tail(5) # You can change to .tail(2) for a quicker peek
  896. # %% [markdown]
  897. # ### ✂️ Data Slicing for G28 Stress
  899. # %%
  900. # Extract signal data (excluding metadata)
  901. cells_traces_fig1_stress_tot = data_fig1_stress.iloc[1:, 1:]
  903. # Get shape info
  904. n_cells_fig1_stress, _ = cells_traces_fig1_stress_tot.shape
  906. # Preview bottom rows
  907. cells_traces_fig1_stress_tot.tail()
  909. # %%
  910. # Convert DataFrame to NumPy array for numerical processing
  911. data_fig1_stress_tot = cells_traces_fig1_stress_tot.values.astype(float)
  913. # Slice the data before stimulation onset
  914. data_fig1_stress_1 = data_fig1_stress_tot[:, :stress_onset]
  916. # Slice data during stimulation period
  917. data_fig1_stress_2 = data_fig1_stress_tot[:, stress_onset:stress_offset]
  919. # %%
  920. # Rescale pre-stimulation traces using min-max normalisation (row-wise)
  921. normalised_data_fig1_stress_1 = minmax_scale(data_fig1_stress_1, axis=1)
  923. # Rescale stimulation-period traces using min-max normalisation (row-wise)
  924. normalised_data_fig1_stress_2 = minmax_scale(data_fig1_stress_2, axis=1)
  926. # %% [markdown]
  927. # ### 📈 corrcoef: Standard Pearson Correlation Coefficient
  929. # %%
  930. # Compute the correlation matrix for pre-stimulation traces
  931. max_corr_fig1_stress_1 = np.corrcoef(normalised_data_fig1_stress_1, rowvar=True)
  933. # Compute the correlation matrix for stimulation-period traces
  934. max_corr_fig1_stress_2 = np.corrcoef(normalised_data_fig1_stress_2, rowvar=True)
  936. # %% [markdown]
  937. # # 🖼️ Figure 1
  939. # %%
  940. # Riemannian distances between empirical correlation matrices and a null model
  941. # - Each value reflects how far a real correlation matrix is from a distribution of random matrices
  942. # - Random matrices generated via Metropolis-Hastings sampling
  943. # - Distances are normalised by the mean pairwise distance among sampled matrices
  945. # Control condition: 12 samples across 4 animals (G33–G36)
  946. control_dist = [0.6484, 0.6310, 0.6835, 0.3532, 0.4503, 0.4678, 0.4226, 0.4540, 0.5091, 0.4670, 0.5697, 0.5820]
  948. control_animal = ['G33', 'G33', 'G33', 'G34', 'G34', 'G34', 'G35', 'G35', 'G35', 'G36', 'G36', 'G36']
  950. # Stimulated condition: 20 samples across 6 animals (G7–G9, G28–G29, G37)
  951. stim_dist = [0.8543, 0.8024, 0.8918, 0.7506, 0.8326, 0.8286, 0.7508, 0.7120, 0.8152, 0.7978,
  952. 0.5823, 0.7944, 0.6782, 0.4742, 0.3286, 0.4824, 0.7652, 0.7087, 0.6085, 0.6523]
  954. stim_animal = ['G7', 'G7', 'G8', 'G8', 'G8', 'G8', 'G9', 'G9', 'G9', 'G9',
  955. 'G28', 'G28','G28', 'G29', 'G29', 'G29', 'G37', 'G37', 'G37', 'G37']
  957. # Stress condition: 7 samples across 6 animals (G7–G9, G28–G29, G37)
  958. # - Note: G37 contributes two samples from distinct sessions
  959. stress_dist = [1.001113156, 1.002419683, 0.969890523, 0.7640226, 0.790428686, 0.885770891, 0.810616632]
  961. stress_animal = ['G7', 'G8', 'G9', 'G28', 'G29', 'G37', 'G37']
  963. # %%
  964. # AUC values represent integrated response over time windows
  965. # - Each list contains per-sample AUCs for a given condition and time window
  966. # - Negative values may reflect suppression or baseline drift
  967. # - Stress values show large positive outliers, suggesting high activation
  969. # ⏱ 30-second window
  970. AUC_stim_30 = [-3.2941, -0.8233, 1.4743, -1.9491, 2.5396, 1.0696, -0.6789, -0.4772, -0.6625, -2.1106,
  971. -2.6996, 10.156, 1.2475, 4.7087, -1.5431, 3.906, 1.4121, -2.2647, -1.1277, -3.9984]
  972. AUC_stress_30 = [3.9668, 4.511, 2.5152, 60.6365, 23.0652, 22.4978, 41.8574]
  974. # ⏱ 60-second window
  975. AUC_stim_60 = [-3.0155, -1.267, 1.6671, -2.3267, 0.952, -0.4225, -2.2116, 1.5772, -0.1854, -1.7976,
  976. -1.8799, 14.8400, -2.0294, 4.1857, -0.9583, 1.1032, 0.0982, -2.7538, -0.5514, -3.5729]
  977. AUC_stress_60 = [4.7454, 2.3182, 0.4937, 50.6340, 21.3518, 26.5827, 51.8101]
  979. # ⏱ 90-second window
  980. AUC_stim_90 = [-1.2230, -1.1646, 1.5550, -1.1195, 1.0739, -0.7815, -2.6634, 1.4759, -0.1702, -3.1705,
  981. -0.7590, 13.0363, -2.0926, 3.4858, 1.1061, 2.0930, -2.0750, -4.4860, -1.3451, -3.1809]
  982. AUC_stress_90 = [4.2257, 1.6342, -1.6458, 39.6093, 17.7396, 24.2487, 50.9404]
  984. # ⏱ 120-second window
  985. AUC_stim_120 = [-0.6409, -0.7567, 0.7624, -0.5384, 0.7329, -0.8672, 0.0113, 1.1853, 0.6831, -0.7784,
  986. -0.8053, 5.1225, -0.5467, 1.7969, -0.0616, 0.3441, -1.2354, -2.4234, -0.7771, -1.499]
  987. AUC_stress_120 = [1.0498, -1.5431, -0.7729, 17.6906, 11.2251, 6.8057, 20.8275]
  989. # %%
  990. # --- Helper function to build sub-DataFrame ---
  991. def build_df(animal_list, values, condition):
  992. df = pd.DataFrame({
  993. 'Animal_id': animal_list,
  994. 'y': values
  995. })
  996. # Add condition and compute repeat count per animal
  997. df['condition'] = condition
  998. df['repeat'] = df.groupby('Animal_id').cumcount() + 1
  999. return df[['Animal_id', 'condition', 'repeat', 'y']]
  1002. # --- Combine all conditions ---
  1003. df_control = build_df(control_animal, control_dist, 'control')
  1004. df_stim = build_df(stim_animal, stim_dist, 'stim')
  1005. df_stress = build_df(stress_animal, stress_dist, 'stress')
  1007. # Merge into one DataFrame
  1008. df = pd.concat([df_control,df_stim, df_stress], ignore_index=True)
  1010. # Sort for readability
  1011. df = df.sort_values(['condition', 'Animal_id', 'repeat']).reset_index(drop=True)
  1013. print(df)
  1015. # %%
  1016. # Nested repeats within animal:
  1017. # y ~ condition + (1 | Animal_id)
  1018. m_nested = smf.mixedlm(
  1019. "y ~ C(condition, Treatment(reference='control'))",
  1020. df,
  1021. groups=df["Animal_id"] # (1 | Animal_id)
  1022. ).fit(reml=True, method="lbfgs", maxiter=600, disp=False)
  1023. print("\n=== Nested model (animal + animal:repeat) ===")
  1024. print(m_nested.summary())
  1025. pvalues_variable = m_nested.pvalues
  1027. # %%
  1028. # Extract p-values for Stim and Stress vs Control
  1029. p_control_Stim = m_nested.pvalues.get("C(condition, Treatment(reference='control'))[T.stim]", None)
  1030. p_control_Stress = m_nested.pvalues.get("C(condition, Treatment(reference='control'))[T.stim]", None)
  1031. print(p_control_Stim, p_control_Stress)
  1033. # %%
  1034. # Compute p-value for stim vs stress using t-test
  1035. # Extract coefficients and standard errors
  1036. params = m_nested.params
  1037. bse = m_nested.bse
  1039. # Coefficients for stim and stress
  1040. b_stim = params["C(condition, Treatment(reference='control'))[T.stim]"]
  1041. b_stress = params["C(condition, Treatment(reference='control'))[T.stress]"]
  1043. # Calculate the difference in coefficients for stress and stim conditions
  1044. diff = b_stim - b_stress
  1046. # Retrieve the covariance matrix of the fixed effects, which contains the variances and covariances of the estimated coefficients.
  1047. cov_matrix = m_nested.cov_params()
  1049. # Standard error of the difference
  1050. se_diff = np.sqrt(cov_matrix.loc["C(condition, Treatment(reference='control'))[T.stim]", "C(condition, Treatment(reference='control'))[T.stim]"] +
  1051. cov_matrix.loc["C(condition, Treatment(reference='control'))[T.stress]", "C(condition, Treatment(reference='control'))[T.stress]"] -
  1052. 2 * cov_matrix.loc["C(condition, Treatment(reference='control'))[T.stim]", "C(condition, Treatment(reference='control'))[T.stress]"])
  1054. # t-statistic
  1055. t_stat = diff / se_diff
  1057. # Degrees of freedom (approximation)
  1058. df = m_nested.df_resid
  1060. # p-value (two-tailed)
  1061. p_Stress_Stim = 2 * t.sf(np.abs(t_stat), df)
  1063. print(p_Stress_Stim)
  1065. # %%
  1066. # Diverging colour palette (brown → teal), normalised to [0, 1] for plotting
  1067. mycolors = np.array([
  1068. [140, 81, 10],
  1069. [216, 179, 101],
  1070. [246, 232, 195],
  1071. [199, 234, 229],
  1072. [ 90, 180, 172],
  1073. [ 1, 102, 94]
  1074. ]) / 255.0
  1076. # Reduced colour palette for control conditions
  1077. colors_ctrl = np.array([
  1078. [166, 97, 26], # #a6611a
  1079. [223, 194, 125], # #dfc27d
  1080. [128, 205, 193], # #80cdc1
  1081. [ 1, 133, 113] # #018571
  1082. ]) / 255.0
  1084. # %%
  1085. # Create a figure with specific dimensions
  1086. fig = plt.figure(figsize=(15, 9), dpi=600)
  1088. # Set up the grid layout for subplots
  1089. gs = gridspec.GridSpec(3, 8, figure=fig)
  1091. # 🔹 Plot 1: Stimulus Data (Baseline)
  1092. ax = fig.add_subplot(gs[0:4])
  1093. # Normalise stimulus data: z-score across rows (cells), then min-max across columns (time)
  1094. stim_data = minmax_scale(data_fig1_stim_tot[:, :stim_offset], axis=1)
  1095. plt.imshow(stim_data, cmap="viridis", interpolation='none')
  1097. # Set axis limits and ticks for clarity
  1098. ax.set_xlim(0, stim_offset)
  1099. ax.set_ylim(0, n_cells_fig1_stim-1)
  1100. ax.set_xticks(np.arange(0, 2500, 500))
  1101. ax.set_xticklabels(np.arange(0, 250, 50))
  1102. ax.set_yticks([_ for _ in list(range(n_cells_fig1_stim)) if _ % 10 == 0])
  1103. ax.set_yticklabels([_ for _ in list(range(n_cells_fig1_stim)) if _ % 10 == 0])
  1104. ax.set_ylabel('cell ID')
  1105. ax.set_xlabel('time [sec]')
  1106. ax.set_aspect('auto')
  1108. # Add colourbar for stimulus data
  1109. cmap = plt.get_cmap("viridis")
  1110. sm = ScalarMappable(cmap=cmap)
  1111. sm.set_clim(vmin=0, vmax=1)
  1112. cb = plt.colorbar(cmap=cmap, fraction=0.04, pad=0.03)
  1113. cb.set_label('$\Delta$ F$/$F$_0$ (scaled)')
  1115. # 🔹 Plot 2: Stress Data (Stress Onset)
  1116. ax = fig.add_subplot(gs[4:8])
  1117. # Normalise stress data: z-score across rows, then min-max across columns
  1118. stress_data = minmax_scale(data_fig1_stress_tot[:, :stress_offset], axis=1)
  1119. plt.imshow(stress_data, cmap="viridis", interpolation='none')
  1121. # Set axis limits and ticks for clarity
  1122. ax.set_xlim(0, stress_offset)
  1123. ax.set_ylim(0, n_cells_fig1_stress-1)
  1124. ax.set_xticks(np.arange(0, 4000, 500))
  1125. ax.set_xticklabels(np.arange(0, 400, 50))
  1126. ax.set_yticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1127. ax.set_yticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1128. ax.set_ylabel('cell ID')
  1129. ax.set_xlabel('time [sec]')
  1130. ax.set_aspect('auto')
  1132. # Add colourbar for stress data
  1133. cmap = plt.get_cmap("viridis")
  1134. sm = ScalarMappable(cmap=cmap)
  1135. sm.set_clim(vmin=0, vmax=1)
  1136. cb = plt.colorbar(cmap=cmap, fraction=0.04, pad=0.03)
  1137. cb.set_label('$\Delta$ F$/$F$_0$ (scaled)')
  1139. # 🔹 Plot 3: Heatmap for Max Correlation (Baseline)
  1140. ax = fig.add_subplot(gs[8:10])
  1141. # Heatmap of pairwise correlations during baseline (stim)
  1142. sns.heatmap(max_corr_fig1_stim_1, cbar_kws={'label': 'Corr'}, cmap="viridis", vmin=-1, vmax=1)
  1143. ax.set_title('Baseline')
  1144. ax.set(xlabel='cell ID', ylabel='cell ID')
  1145. ax.set_xticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1146. ax.set_xticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1147. ax.set_yticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1148. ax.set_yticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1149. ax.invert_yaxis()
  1151. # Add outline around the whole plot
  1152. for spine in ax.spines.values():
  1153. spine.set_visible(True)
  1154. spine.set_linewidth(1.5)
  1155. spine.set_edgecolor('gray')
  1157. # Add grey outline around the colourbar
  1158. cbar = ax.collections[0].colorbar
  1159. cbar.outline.set_visible(True)
  1160. cbar.outline.set_edgecolor('grey')
  1161. cbar.outline.set_linewidth(1.0)
  1163. # 🔹 Plot 4: Heatmap for Max Correlation (Stimulation)
  1164. ax = fig.add_subplot(gs[10:12])
  1165. # Heatmap of pairwise correlations during stimulation
  1166. sns.heatmap(max_corr_fig1_stim_2, cbar_kws={'label': 'Corr'}, cmap="viridis", vmin=-1, vmax=1)
  1167. ax.set_title('Stimulation')
  1168. ax.set(xlabel='cell ID')
  1169. ax.set_xticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1170. ax.set_xticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1171. ax.set_yticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1172. ax.set_yticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1173. ax.invert_yaxis()
  1175. # Add outline around the whole plot
  1176. for spine in ax.spines.values():
  1177. spine.set_visible(True)
  1178. spine.set_linewidth(1.5)
  1179. spine.set_edgecolor('gray')
  1181. # Add grey outline around the colourbar
  1182. cbar = ax.collections[0].colorbar
  1183. cbar.outline.set_visible(True)
  1184. cbar.outline.set_edgecolor('grey')
  1185. cbar.outline.set_linewidth(1.0)
  1187. # 🔹 Plot 5: Heatmap for Max Correlation (Stress Baseline)
  1188. ax = fig.add_subplot(gs[12:14])
  1189. # Heatmap of pairwise correlations during stress baseline
  1190. sns.heatmap(max_corr_fig1_stress_1, cbar_kws={'label': 'Corr'}, cmap="viridis", vmin=-1, vmax=1)
  1191. ax.set_title('Baseline')
  1192. ax.set(xlabel='cell ID', ylabel='cell ID')
  1193. ax.set_xticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1194. ax.set_xticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1195. ax.set_yticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1196. ax.set_yticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1197. ax.invert_yaxis()
  1199. # Add outline around the whole plot
  1200. for spine in ax.spines.values():
  1201. spine.set_visible(True)
  1202. spine.set_linewidth(1.5)
  1203. spine.set_edgecolor('gray')
  1205. # Add grey outline around the colourbar
  1206. cbar = ax.collections[0].colorbar
  1207. cbar.outline.set_visible(True)
  1208. cbar.outline.set_edgecolor('grey')
  1209. cbar.outline.set_linewidth(1.0)
  1211. # 🔹 Plot 6: Heatmap for Max Correlation (Stress)
  1212. ax = fig.add_subplot(gs[14:16])
  1213. # Heatmap of pairwise correlations during stress period
  1214. sns.heatmap(max_corr_fig1_stress_2, cbar_kws={'label': 'Corr'}, cmap="viridis", vmin=-1, vmax=1)
  1215. ax.set_title('Stress')
  1216. ax.set(xlabel='cell ID')
  1217. ax.set_xticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1218. ax.set_xticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1219. ax.set_yticks([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1220. ax.set_yticklabels([_ for _ in list(range(n_cells_fig1_stress)) if _ % 10 == 0])
  1221. ax.invert_yaxis()
  1223. # Add outline around the whole plot
  1224. for spine in ax.spines.values():
  1225. spine.set_visible(True)
  1226. spine.set_linewidth(1.5)
  1227. spine.set_edgecolor('gray')
  1229. # Add grey outline around the colourbar
  1230. cbar = ax.collections[0].colorbar
  1231. cbar.outline.set_visible(True)
  1232. cbar.outline.set_edgecolor('grey')
  1233. cbar.outline.set_linewidth(1.0)
  1235. # 🔹 Plot 7: Swarm plot for stimulation p-values
  1236. colour_list = ['blue', 'red', 'green', 'orange', 'purple', 'brown']
  1237. ax = fig.add_subplot(gs[16:19])
  1238. # Violin + swarm plot for control condition
  1239. cc1 = control_dist
  1241. control_animals = pd.unique(control_animal)
  1242. control_palette = dict(zip(
  1243. control_animals,
  1244. colors_ctrl[:len(control_animals)]
  1245. ))
  1247. df = pd.DataFrame(dict(x=np.repeat([0], len(cc1)), y=cc1, c=control_animal))
  1248. sns.violinplot(x="x", y="y", data=df, order=['0','1','2'], fill=False, color='black')
  1249. sns.swarmplot(
  1250. x="x",
  1251. y="y",
  1252. hue="c",
  1253. data=df,
  1254. palette=control_palette,
  1255. size=7,
  1256. marker='s',
  1257. legend=False
  1258. )
  1260. # Violin + swarm plot for stimulation condition
  1262. animals = pd.unique(stim_animal)
  1263. palette = dict(zip(animals, mycolors))
  1265. cc2 = stim_dist
  1267. df = pd.DataFrame(dict(x=np.repeat([1], len(cc2)), y=cc2, c=stim_animal))
  1268. sns.violinplot(x="x", y="y", data=df, order=['0','1','2'], fill=False, color='black')
  1269. sns.swarmplot(
  1270. x="x",
  1271. y="y",
  1272. hue="c",
  1273. data=df,
  1274. palette=palette,
  1275. order=[1], # x is numeric and always 1
  1276. size=7,
  1277. legend=False
  1278. )
  1280. # Violin + swarm plot for stress condition
  1281. cc3 = stress_dist
  1282. df = pd.DataFrame(dict(
  1283. x=np.repeat(2, len(cc3)),
  1284. y=cc3,
  1285. c=stress_animal
  1286. ))
  1287. sns.violinplot(x="x", y="y", data=df, order=['0','1','2'], fill=False, color='black')
  1288. sns.swarmplot(
  1289. x="x",
  1290. y="y",
  1291. hue="c",
  1292. data=df,
  1293. palette=palette,
  1294. order=[2],
  1295. size=7,
  1296. legend=False
  1297. )
  1299. # Add statistical annotations between groups
  1300. # Control vs Stim
  1301. ax.plot([0.1, 0.9], [1.2, 1.2], 'k-', lw=1)
  1302. ax.plot([0.1, 0.1], [1.2, 1.18], 'k-', lw=1)
  1303. ax.plot([0.9, 0.9], [1.2, 1.18], 'k-', lw=1)
  1304. _, kruskal_p_value = kruskal(control_dist, stim_dist)
  1305. ax.annotate(pval_to_star(p_control_Stim), (1/2, 1.21), ha='center')
  1307. # Stim vs Stress
  1308. ax.plot([1.1, 1.9], [1.2, 1.2], 'k-', lw=1)
  1309. ax.plot([1.1, 1.1], [1.2, 1.18], 'k-', lw=1)
  1310. ax.plot([1.9, 1.9], [1.2, 1.18], 'k-', lw=1)
  1311. _, kruskal_p_value = kruskal(stress_dist, stim_dist)
  1312. ax.annotate(pval_to_star(p_Stress_Stim), (3/2, 1.21), ha='center')
  1314. # Control vs Stress
  1315. ax.plot([0.1, 1.9], [1.35, 1.35], 'k-', lw=1)
  1316. ax.plot([0.1, 0.1], [1.35, 1.33], 'k-', lw=1)
  1317. ax.plot([1.9, 1.9], [1.35, 1.33], 'k-', lw=1)
  1318. _, kruskal_p_value = kruskal(control_dist, stress_dist)
  1319. ax.annotate(pval_to_star(p_control_Stress), (1, 1.36), ha='center')
  1321. # Final axis settings
  1322. ax.set_ylim(0.0, 1.5)
  1323. ax.set_xlabel('')
  1324. ax.set_ylabel('Distance')
  1325. ax.set_xticks(range(3), labels=['Control', 'Stimulation', 'Stress'])
  1327. # 🔹 Plot 8: Swarm plot for stress p-values
  1328. ax = fig.add_subplot(gs[19:23])
  1330. # Time windows used for AUC calculation
  1331. time_points = [30, 60, 90, 120]
  1333. # Grouped AUC data for stim and stress conditions
  1334. stim_data = [AUC_stim_30, AUC_stim_60, AUC_stim_90, AUC_stim_120]
  1335. stress_data = [AUC_stress_30, AUC_stress_60, AUC_stress_90, AUC_stress_120]
  1337. # Compute median and SEM for each time window
  1338. stim_medians = [np.median(d) for d in stim_data]
  1339. stim_sem = [np.std(d, ddof=1) / np.sqrt(len(d)) for d in stim_data]
  1341. stress_medians = [np.median(d) for d in stress_data]
  1342. stress_sem = [np.std(d, ddof=1) / np.sqrt(len(d)) for d in stress_data]
  1344. # Create a DataFrame for plotting
  1345. df = pd.DataFrame({
  1346. 'Time': time_points,
  1347. 'Stim_Median': stim_medians,
  1348. 'Stim_SEM': stim_sem,
  1349. 'Stress_Median': stress_medians,
  1350. 'Stress_SEM': stress_sem
  1351. })
  1355. ## Create the plot for AUC data
  1357. # -----------------------------
  1358. # Positions
  1359. # -----------------------------
  1360. x = np.arange(len(time_points))
  1361. offset = 0.18
  1363. stim_pos = x - offset
  1364. stress_pos = x + offset
  1366. # -----------------------------
  1367. # Plot
  1368. # -----------------------------
  1371. bp_stim = plt.boxplot(
  1372. stim_data,
  1373. positions=stim_pos,
  1374. widths=0.3,
  1375. patch_artist=True,
  1376. medianprops=dict(color='tab:blue', linewidth=2),
  1377. flierprops=dict(
  1378. marker='o',
  1379. markersize=5,
  1380. markerfacecolor='tab:blue',
  1381. markeredgecolor='tab:blue',
  1382. alpha=0.6
  1383. )
  1384. )
  1386. bp_stress = plt.boxplot(
  1387. stress_data,
  1388. positions=stress_pos,
  1389. widths=0.3,
  1390. patch_artist=True,
  1391. medianprops=dict(color='tab:orange', linewidth=2),
  1392. flierprops=dict(
  1393. marker='o',
  1394. markersize=5,
  1395. markerfacecolor='tab:orange',
  1396. markeredgecolor='tab:orange',
  1397. alpha=0.6
  1398. )
  1399. )
  1401. # -----------------------------
  1402. # Box transparency & colors
  1403. # -----------------------------
  1404. for box in bp_stim['boxes']:
  1405. box.set(facecolor='tab:blue', alpha=0.45)
  1407. for box in bp_stress['boxes']:
  1408. box.set(facecolor='tab:orange', alpha=0.45)
  1410. # -----------------------------
  1411. # Connect medians
  1412. # -----------------------------
  1413. stim_medians = [np.median(d) for d in stim_data]
  1414. stress_medians = [np.median(d) for d in stress_data]
  1416. ax.plot(stim_pos, stim_medians, 'o', color='tab:blue', linewidth=2, label='Stim')
  1417. ax.plot(stress_pos, stress_medians, 'o', color='tab:orange', linewidth=2, label='Stress')
  1422. # Plot median ± SEM for stimulation
  1423. #ax.errorbar(df['Time'], df['Stim_Median'], yerr=df['Stim_SEM'],
  1424. # label='Stim', marker='o', capsize=5, linestyle='-')
  1426. # Plot median ± SEM for stress
  1427. #ax.errorbar(df['Time'], df['Stress_Median'], yerr=df['Stress_SEM'],
  1428. # label='Stress', marker='s', capsize=5, linestyle='--')
  1430. # Annotate significance for each time window using permulation test (see data/calcium_imaging/figures/fig1_stim_stress_permutation_test.py)
  1431. mw_p_value = 0.0132
  1432. ax.annotate(pval_to_star(mw_p_value), (0, 70), ha='center')
  1434. mw_p_value = 0.0121
  1435. ax.annotate(pval_to_star(mw_p_value), (1, 70), ha='center')
  1437. mw_p_value = 0.0421
  1438. ax.annotate(pval_to_star(mw_p_value), (2, 70), ha='center')
  1440. mw_p_value = 0.09709
  1441. ax.annotate(pval_to_star(mw_p_value), (3, 70), ha='center')
  1443. # Final axis settings
  1444. ax.set_ylim(-10, 80)
  1445. ax.set_xlabel('Time [sec]')
  1446. ax.set_ylabel('Time-Averaged\n Activity')
  1447. ax.set_xticks([0,1,2,3], labels=time_points)
  1448. ax.legend(
  1449. )
  1453. # Adjust layout to prevent overlap and ensure clean spacing
  1454. plt.tight_layout()
  1457. # Display the figure
  1458. # fig.savefig("Fig1_revised.svg", format="svg", bbox_inches="tight", dpi=600)
  1459. plt.show()
  1461. # %% [markdown]
  1462. # # 🧪 Full Workflow: G28 Stimulation (Figure 2)
  1464. # %% [markdown]
  1465. # An example of data analysis pipline for stimulation experimental setup.
  1467. # %%
  1468. # Load and transpose the data
  1469. data_fig2_stim = pd.read_csv(stim_fig2_csv_fname).T
  1471. # Clean up row labels
  1472. data_fig2_stim.index = data_fig2_stim.index.str.strip()
  1474. # Preview the last 6 cells
  1475. data_fig2_stim.tail(6)
  1477. # %% [markdown]
  1478. # ## ✂️ Data Slicing for G28 Stimulation
  1480. # %%
  1481. # Extract signal data (excluding metadata)
  1482. cells_traces_tot_old = data_fig2_stim.iloc[1:, 1:]
  1484. # Get shape info
  1485. total_cells_stim, data_stim_length_tot = cells_traces_tot_old.shape
  1487. # Preview bottom rows
  1488. cells_traces_tot_old.tail()
  1490. # %%
  1491. # Extract time axis from first row (excluding first column)
  1492. ts_stim_tot = np.asarray(data_fig2_stim.iloc[0, 1:], dtype=float)
  1494. # Time between frames (first timestamp is 0)
  1495. sec_per_frame = ts_stim_tot[1]
  1497. # Sampling frequency (Hz)
  1498. fs = 1 / sec_per_frame
  1500. # %%
  1501. # Convert fluorescence traces to NumPy array
  1502. data_stim_artifact = np.asarray(cells_traces_tot_old, dtype=float)
  1504. # %% [markdown]
  1505. # ## 🔍 Finding UCN3 Cells
  1507. # %% [markdown]
  1508. # ### 🌀 1. Power Spectral Density (PSD) Filtering
  1510. # %%
  1511. # Define the stimulus frequency and initialise a list to store labels of removed cells based on PSD threshold
  1512. stimulus_freq = 0.1 # Stimulus frequency in Hz
  1513. labels_removed_PSD = [] # List to store indices of cells with high PSD at stimulus frequency
  1515. # Loop through each cell trace to assess frequency-domain response
  1516. for i in range(total_cells_stim):
  1517. # Extract stimulation-period trace for cell i
  1518. trace_segment = data_stim_artifact[i, stim_onset:stim_offset]
  1520. # Compute power spectral density using periodogram
  1521. freq, Pxx_den = periodogram(trace_segment, fs=fs)
  1523. # Normalise PSD to its peak value for comparability
  1524. Pxx_den /= np.max(Pxx_den)
  1526. # Identify index corresponding to stimulus frequency
  1527. stimulus_freq_idx = round(stimulus_freq * (stim_offset - stim_onset) / fs)
  1529. # Thresholding: retain cells with strong PSD at stimulus frequency
  1530. if Pxx_den[stimulus_freq_idx] > PSD_thr:
  1531. labels_removed_PSD.append(i)
  1533. # Output: indices of cells flagged for removal based on PSD criterion
  1534. print("Indices of cells removed based on PSD threshold:", labels_removed_PSD)
  1536. # %%
  1537. # Create binary mask to flag artifact-affected cells
  1538. # - 0: cell identified as artifact (high PSD at stimulus frequency)
  1539. # - 1: cell retained for further analysis
  1540. labels_artifact = np.asarray([
  1541. 0 if cell in labels_removed_PSD else 1
  1542. for cell in range(total_cells_stim)
  1543. ])
  1545. # Optional: Print the binary mask to verify which cells were flagged
  1546. # 0 = artifact (failed PSD threshold), 1 = valid (passed PSD threshold)
  1547. print(labels_artifact)
  1549. # %% [markdown]
  1550. # ### 🔢 2. K-Means Clustering
  1552. # %%
  1553. # Apply scikit-learn K-means clustering to sort activity traces
  1555. # Step 1: Normalise the activity data across time (columns) to ensure fair clustering
  1556. scaled_data = scale(data_stim_artifact[:, stim_onset:stim_offset], axis=0)
  1558. # Step 2: Set the number of clusters (2 groups: responsive vs non-responsive)
  1559. Ks = 2
  1561. # Step 3: Initialise and run K-means clustering
  1562. # - random_state ensures reproducibility
  1563. # - n_init = 10 means the algorithm runs 10 times with different centroid seeds
  1564. # - max_iter = 100 limits the number of iterations per run
  1565. kmeans = KMeans(n_clusters=Ks, random_state=11111, n_init=10, max_iter=100)
  1567. # Step 4: Fit the model to the scaled data and assign cluster labels to each cell
  1568. labels_artifact_KM = kmeans.fit_predict(scaled_data)
  1570. # Step 5: Sort cell indices based on their assigned cluster labels
  1571. # This helps organise or visualise cells by their activity patterns
  1572. i_labels_artifact = np.argsort(labels_artifact_KM)
  1574. # Step 6: Display the sorted indices for inspection or downstream analysis
  1575. print(f"Sorted indices of cells by cluster ID:\n{i_labels_artifact}")
  1577. # %%
  1578. # Determine which cluster to remove based on relative size
  1580. # Step 1: Count the number of cells in each cluster
  1581. size_cluster_0 = np.sum(labels_artifact_KM == 0)
  1582. size_cluster_1 = np.sum(labels_artifact_KM == 1)
  1584. # Step 2: Compare cluster sizes to identify the smaller one
  1585. # The assumption here is that the smaller cluster likely represents artifacts or noise
  1586. if size_cluster_0 < size_cluster_1:
  1587. cluster_removed = 0 # Smaller cluster is flagged for removal
  1588. cluster_kept = 1 # Larger cluster is retained
  1589. print("Cluster #0 (smaller) is removed.")
  1590. elif size_cluster_1 < size_cluster_0:
  1591. cluster_removed = 1
  1592. cluster_kept = 0
  1593. print("Cluster #1 (smaller) is removed.")
  1594. else:
  1595. # Edge case: one cluster has zero members, so no valid removal decision
  1596. print("Clusters are equal in size — no removal applied.")
  1598. # %%
  1599. # Update artifact labels based on K-means clustering results
  1601. # Step 1: Check if the cluster marked for removal contains any cells
  1602. if np.any(labels_artifact_KM == cluster_removed):
  1603. # Step 2: Iterate through all cells
  1604. for i in range(len(data_stim_artifact)):
  1605. # Only update cells that were previously marked as valid (1)
  1606. # and now belong to the cluster flagged for removal
  1607. if labels_artifact[i] == 1 and labels_artifact_KM[i] == cluster_removed:
  1608. labels_artifact[i] = 0 # Mark as artifact
  1610. # Step 3: Explicitly define which cluster was removed and which was kept
  1611. cluster_removed = 0
  1612. cluster_kept = 1
  1614. # Step 4: Extract indices of cells now marked as artifacts
  1615. labels_removed = np.where(labels_artifact == cluster_removed)[0]
  1617. # Optional: print summary for verification
  1618. print(f"Number of (UCN) cells removed: {len(labels_removed)}")
  1619. print(f"Indices of removed cells: {labels_removed}")
  1621. # %%
  1622. # Remove resonated (artifact) cells from fluorescence intensity traces
  1624. # Step 1: Get the row labels (e.g., cell IDs) corresponding to removed indices
  1625. removed_row_labels = cells_traces_tot_old.index[labels_removed]
  1627. # Step 2: Drop the identified artifact cells from the original DataFrame
  1628. # - axis=0 means rows are being removed
  1629. # - inplace=False ensures the original DataFrame remains unchanged
  1630. cells_traces_stim_tot = cells_traces_tot_old.drop(
  1631. index=removed_row_labels, axis=0, inplace=False
  1632. )
  1634. # Step 3: Count the number of remaining (valid) cells after removal
  1635. n_cells_stim = len(cells_traces_stim_tot)
  1637. # Optional: print summary
  1638. print(f"Number of GABA cells retained after artifact removal: {n_cells_stim}")
  1640. # %% [markdown]
  1641. # ## Minmax Scaling over each trace
  1643. # %%
  1644. # Convert fluorescence traces to NumPy array (float)
  1645. data_stim_tot = np.asarray(cells_traces_stim_tot, dtype=float)
  1647. # Apply Min-Max scaling to each trace individually
  1648. # - Rescales each cell's activity to [0, 1] across its full time course
  1649. normalised_data_stim_tot = minmax_scale(data_stim_tot, axis=1)
  1651. # %% [markdown]
  1652. # ## ⚡️ **G28 Stimulation**
  1653. #
  1654. # Hereafter, I only focus on the data while the stimulation is applied at 2 minutes after starting recording with 0.1 Hz (5 s on followed by 5 s off).
  1656. # %% [markdown]
  1657. # ### ✂️ Data Slicing
  1659. # %%
  1660. # Slice fluorescence data to isolate stimulation period
  1661. cells_traces_stim_2 = cells_traces_stim_tot.iloc[:, stim_onset:stim_offset]
  1663. # Min-max normalised fluorescence during stimulation
  1664. normalised_data_stim_2 = normalised_data_stim_tot[:, stim_onset:stim_offset]
  1666. # Update stimulation-period length
  1667. data_stim_length_2 = cells_traces_stim_2.shape[1]
  1669. # Optional: preview first few traces during stimulation
  1670. cells_traces_stim_2.head()
  1672. # %%
  1673. # Create time axis for stimulation period
  1674. # - Ensures alignment with sliced fluorescence data
  1675. ts_stim_2 = ts_stim_tot[stim_onset:stim_offset]
  1677. # %%
  1678. # Convert stimulation-period fluorescence data to NumPy array (float)
  1679. # - Enables numerical operations like clustering, filtering, or plotting
  1680. data_stim_2 = np.asarray(cells_traces_stim_2, dtype=float)
  1682. # %% [markdown]
  1683. # ### Clustering
  1685. # %% [markdown]
  1686. # #### 🌲 Hierarchical (corrcoef)
  1688. # %%
  1689. # Compute pairwise correlation-based distances
  1690. # - Distance metric: 1 - Pearson correlation coefficient
  1691. # - pdist returns a condensed distance matrix (n*(n-1)/2 elements)
  1692. dist_corr_stim_2 = pdist(normalised_data_stim_2,
  1693. lambda x, y: 1 - np.corrcoef(x, y)[0, 1])
  1695. # %%
  1696. # Convert condensed distance matrix to square form
  1697. # - Required for AgglomerativeClustering with 'precomputed' metric
  1698. dist_matrix = squareform(dist_corr_stim_2)
  1700. # Apply hierarchical clustering to sort activity traces
  1701. # - AgglomerativeClustering performs hierarchical clustering.
  1702. # - `n_clusters=Ks_2` specifies the number of clusters.
  1703. # - `metric='precomputed'` indicates that the distance matrix has already been computed.
  1704. # - `linkage='average'` uses the average linkage criterion for merging clusters.
  1705. # - `compute_distances=True` ensures that pairwise distances are retained.
  1706. ac = AgglomerativeClustering(
  1707. n_clusters=Ks_2,
  1708. metric='precomputed',
  1709. linkage='average',
  1710. compute_distances=True
  1711. ).fit(dist_matrix)
  1713. # Retrieve the cluster labels assigned to each trace
  1714. labels_hier_stim_2 = ac.labels_
  1716. # %%
  1717. # Compute pairwise Pearson correlation matrix for stimulation data
  1718. # - Avoids redundant diagonal comparisons
  1719. # - Uses z-scored traces for consistency
  1720. max_corr_stim_2 = np.zeros((n_cells_stim, n_cells_stim))
  1722. for i in range(n_cells_stim):
  1723. for j in range(i + 1, n_cells_stim): # Upper triangle only
  1724. trace_i = scale(normalised_data_stim_2[i, :])
  1725. trace_j = scale(normalised_data_stim_2[j, :])
  1726. corr = np.corrcoef(trace_i, trace_j)[0, 1]
  1727. max_corr_stim_2[i, j] = corr
  1728. max_corr_stim_2[j, i] = corr # Symmetric assignment
  1730. # %%
  1731. from sklearn.metrics import silhouette_score
  1733. # Evaluate clustering quality using silhouette score
  1734. # - Measures how well each trace fits within its assigned cluster
  1735. # - Higher scores indicate better-defined clusters
  1736. sil_score = silhouette_score(
  1737. squareform(dist_corr_stim_2), # Full distance matrix
  1738. labels_hier_stim_2, # Cluster labels
  1739. metric="precomputed" # Use provided distances
  1740. )
  1742. # Display the silhouette score
  1743. print(f'Silhouette score: {sil_score:.3f}')
  1745. # %% [markdown]
  1746. # #### 🔢 K-means (corrcoef)
  1748. # %%
  1749. # Apply scikit-learn K-means clustering to the SLxCorr distance matrix
  1750. # Fix the random seed for reproducibility and set the number of clusters (Ks_2)
  1751. km = KMeans(n_clusters=Ks_2, random_state=11112, n_init=10, max_iter=100)
  1753. # Fit the model on the squared form of the SLxCorr distance matrix
  1754. km.fit(squareform(dist_corr_stim_2))
  1756. # Get the labels for the clusters assigned to each cell
  1757. labels_corr_stim_2 = km.labels_
  1759. # %%
  1760. # Plotting the unsorted and sorted (SLxCorr, K-means) distance matrices
  1761. plt.figure(figsize=(16, 6))
  1763. # Subplot 1: Unsorted SLxCorr distance matrix
  1764. plt.subplot(1, 2, 1)
  1765. sns.heatmap(2 * minmax_scale(squareform(dist_corr_stim_2),axis=1), cmap='BrBG')
  1766. plt.title('Unsorted SLxCorr Distance Matrix')
  1768. # Subplot 2: Sorted SLxCorr distance matrix based on K-means labels
  1769. plt.subplot(1, 2, 2)
  1770. # Sorting the SLxCorr distance matrix based on the K-means cluster labels
  1771. idx = np.argsort(labels_corr_stim_2) # Sort indices by K-means labels
  1772. sorted_cov_dist = squareform(dist_corr_stim_2)[idx, :][:, idx] # Apply sorting
  1773. ax = sns.heatmap(2 * minmax_scale(sorted_cov_dist,axis=1),
  1774. cmap='BrBG', cbar_kws={'label': 'Dissimilarity'})
  1775. ax.set(xlabel='Cell ID', ylabel='Cell ID')
  1776. plt.title('Sorted SLxCorr Distance Matrix (K-means)')
  1778. # Show the plot
  1779. plt.show()
  1781. # %% [markdown]
  1782. # # 🧪 Full Workflow: G9 Stress (Figure 2)
  1784. # %% [markdown]
  1785. # An example of data analysis pipline for restraint stress experimental setup.
  1787. # %%
  1788. # Load stress dataset from CSV and transpose it
  1789. # - Transposition makes each row represent a cell (e.g., for time-series or feature vectors)
  1790. data_fig2_stress = pd.read_csv(stress_fig2_csv_fname).T
  1792. # Clean up cell ID labels by stripping whitespace
  1793. # - Ensures consistent indexing and avoids matching issues
  1794. data_fig2_stress.index = data_fig2_stress.index.str.strip()
  1796. # Display the last few rows of the transposed dataset
  1797. # - Useful for inspecting cell-level data or verifying formatting
  1798. data_fig2_stress.tail()
  1800. # %% [markdown]
  1801. # ## ✂️ Data Slicing
  1803. # %%
  1804. # Extract fluorescence intensity traces from deconvoluted stress data
  1805. # - Skip first row (likely metadata or time vector)
  1806. # - Skip first column (likely cell labels or non-numeric info)
  1807. cells_traces_stress_tot = data_fig2_stress.iloc[1:, 1:]
  1809. # Determine dataset dimensions
  1810. # - n_cells_stress: number of cells (rows)
  1811. # - data_stress_length_tot: number of time points or features (columns)
  1812. n_cells_stress, data_stress_length_tot = cells_traces_stress_tot.shape
  1814. # Preview the last few rows of the fluorescence intensity matrix
  1815. # - Each row = one cell's trace
  1816. # - Each column = one time point or feature
  1817. cells_traces_stress_tot.tail()
  1819. # %%
  1820. # Extract time axis from the first row (excluding first column)
  1821. # - Assumes first row contains timestamps
  1822. # - Converts to NumPy array of floats for numerical operations
  1823. ts_stress_tot = np.asarray(data_fig2_stress.iloc[0, 1:], dtype=float)
  1825. # Determine time per frame (seconds)
  1826. # - Assumes uniform sampling; uses second timestamp as Δt
  1827. sec_per_frame = ts_stress_tot[1]
  1829. # Compute sampling frequency (Hz)
  1830. # - fs = frames per second = 1 / time per frame
  1831. fs = 1 / sec_per_frame
  1833. # %% [markdown]
  1834. # ## Minmax Scaling over each trace
  1836. # %%
  1837. # Convert fluorescence intensity data to NumPy array (float)
  1838. # - Enables fast numerical operations and compatibility with signal processing tools
  1839. # - Each row = one cell's activity trace
  1840. # - Each column = one time point
  1841. data_stress_tot = np.asarray(cells_traces_stress_tot, dtype=float)
  1843. # Normalise each cell's fluorescence trace to [0, 1]
  1844. # - axis=1 ensures scaling is done per row (i.e., per cell)
  1845. # - Preserves temporal dynamics while removing amplitude bias
  1846. normalised_data_stress_tot = minmax_scale(data_stress_tot, axis=1)
  1848. # %% [markdown]
  1849. # ## 🔥 **G7 Stress**
  1850. #
  1851. # Hereafter, I only focus on the data while the stimulation is applied at 2 minutes after starting recording with 0.1 Hz (5 s on followed by 5 s off).
  1853. # %% [markdown]
  1854. # ### ✂️ Data Slicing
  1856. # %%
  1857. # Slice fluorescence intensity traces during stimulation period
  1858. # - Columns correspond to time points; slicing isolates stress window
  1859. cells_traces_stress_2 = cells_traces_stress_tot.iloc[:, stress_onset:stress_offset]
  1861. # Slice normalized traces for the same stimulation window
  1862. # - Ensures amplitude-independent comparison during stress
  1863. normalised_data_stress_2 = normalised_data_stress_tot[:, stress_onset:stress_offset]
  1865. # Get number of time points during stimulation
  1866. data_stress_length_2 = cells_traces_stress_2.shape[1]
  1868. # Preview top rows of sliced raw intensity data
  1869. cells_traces_stress_2.head()
  1871. # %%
  1872. # Extract time axis corresponding to stimulation period
  1873. ts_stress_2 = ts_stress_tot[stress_onset:stress_offset]
  1875. # %%
  1876. # Convert sliced fluorescence intensity data to NumPy array (float)
  1877. data_stress_2 = np.asarray(cells_traces_stress_2.iloc[:, :], dtype=float)
  1879. # %% [markdown]
  1880. # ### Clustering
  1882. # %% [markdown]
  1883. # #### 🌲 Hierarchical (corrcoef)
  1885. # %%
  1886. # Compute pairwise Pearson correlation-based distances
  1887. # - Each trace is standardized (zero mean, unit variance) across time
  1888. # - Distance = 1 - Pearson correlation coefficient
  1889. # → High correlation → small distance
  1890. dist_corr_stress_2 = pdist(
  1891. scale(normalised_data_stress_2, axis=1),
  1892. lambda x, y: 1 - np.corrcoef(x, y)[0, 1]
  1893. )
  1895. # Output is a condensed distance matrix
  1896. # - Suitable for clustering, silhouette scoring, or dendrograms
  1898. # %%
  1899. # Apply hierarchical clustering to standardized activity traces
  1900. # - Uses precomputed Pearson correlation-based distances
  1901. # - Average linkage merges clusters based on mean pairwise distance
  1902. ac = AgglomerativeClustering(
  1903. n_clusters=Ks_2, # 🔢 Desired number of clusters
  1904. metric='precomputed', # 📏 Use custom distance matrix
  1905. linkage='average', # 🔗 Average linkage for merging
  1906. compute_distances=True # 📐 Retain distance info for dendrograms
  1907. ).fit(squareform(dist_corr_stress_2)) # 🔄 Convert condensed to square distance matrix
  1909. # Extract cluster labels for each cell trace
  1910. labels_hier_stress_2 = ac.labels_
  1912. # ✅ `labels_hier_stress_2` contains the cluster assignment for each trace
  1914. # %%
  1915. # Initialize symmetric matrix to store pairwise Pearson correlations
  1916. # - Shape: (n_cells_stress, n_cells_stress)
  1917. # - Diagonal remains zero (self-comparisons skipped)
  1918. max_corr_stress_2 = np.zeros((n_cells_stress, n_cells_stress))
  1920. # Compute pairwise correlations (upper triangle only)
  1921. # - Each trace is standardised (zero mean, unit variance)
  1922. # - Matrix is filled symmetrically to avoid redundant computation
  1923. for first in range(n_cells_stress):
  1924. for second in range(first + 1, n_cells_stress):
  1925. corr = np.corrcoef(
  1926. scale(normalised_data_stress_2[first, :]),
  1927. scale(normalised_data_stress_2[second, :])
  1928. )[0, 1]
  1929. max_corr_stress_2[first, second] = corr
  1930. max_corr_stress_2[second, first] = corr # 🔁 Symmetric fill
  1932. # %%
  1933. # Evaluate clustering quality using silhouette score
  1934. # - Measures how similar each cell is to its own cluster vs other clusters
  1935. # - Uses precomputed Pearson correlation-based distances
  1936. sil_score = silhouette_score(
  1937. squareform(dist_corr_stress_2), # 🔄 Full pairwise distance matrix
  1938. labels_hier_stress_2, # 🏷️ Cluster labels
  1939. metric="precomputed" # 📏 Use custom distance metric
  1940. )
  1942. # Display silhouette score (range: -1 to 1)
  1943. # - Higher = better-defined clusters
  1944. print(f'Silhouette score: {sil_score:.3f}')
  1946. # %% [markdown]
  1947. # #### 🔢 K-means (corrcoef)
  1949. # %%
  1950. # Note: KMeans typically expects feature vectors, not distance matrices.
  1951. # - Applying it to a distance matrix (especially non-Euclidean) may distort clustering.
  1952. # - Consider using MDS or t-SNE to embed the distance matrix into feature space first.
  1954. # Apply K-means clustering to the SLxCorr distance matrix
  1955. # - `random_state` ensures reproducibility
  1956. # - `n_init=10` runs the algorithm multiple times to avoid poor local minima
  1957. # - `max_iter=100` limits iterations per run
  1958. km = KMeans(n_clusters=Ks_2, random_state=11112, n_init=10, max_iter=100)
  1960. # Fit the model on the square distance matrix
  1961. # - Assumes rows are interpretable as feature vectors (approximate)
  1962. km.fit(squareform(dist_corr_stress_2))
  1964. # Extract cluster labels for each cell
  1965. labels_corr_stress_2 = km.labels_
  1967. # %% [markdown]
  1968. # # 🖼️ Figure 2
  1970. # %%
  1971. # Load hierarchical clustering cell percentage data
  1972. # - Each row = one trial
  1973. # - Each column = one cluster category (e.g., Cluster 0, Cluster 1, ...)
  1974. cell_percentage_hier = pd.read_csv(hier_clustering_cell_perc_fname)
  1976. # Preview the bottom rows to inspect structure and completeness
  1977. cell_percentage_hier.tail()
  1979. # %%
  1980. # Extract cluster-1 cell percentages for each stimulated animal
  1981. # NaN values (empty clusters) are removed
  1982. cluster1_percentage = [
  1983. cell_percentage_hier[k].dropna().values
  1984. for k in ['G7 (#1)', 'G8 (#1)', 'G9 (#1)', 'G28 (#1)', 'G29 (#1)', 'G37 (#1)']
  1985. ]
  1987. # Extract cluster-2 cell percentages for each stimulated animal
  1988. # NaN values (empty clusters) are removed
  1989. cluster2_percentage = [
  1990. cell_percentage_hier[k].dropna().values
  1991. for k in ['G7 (#2)', 'G8 (#2)', 'G9 (#2)', 'G28 (#2)', 'G29 (#2)', 'G37 (#2)']
  1992. ]
  1994. flat_cluster1stim = np.concatenate(cluster1_percentage)
  1995. flat_cluster2stim = np.concatenate(cluster2_percentage)
  1997. # Extract cluster percentages for each restraint stress animal
  1998. flat_cluster1stress = cell_percentage_hier['Stress (#1)'].dropna().values
  1999. flat_cluster2stress = cell_percentage_hier['Stress (#2)'].dropna().values
  2001. # Animal labels corresponding to the flattened cluster entries
  2002. stim_animal = [
  2003. 'G7', 'G7',
  2004. 'G8', 'G8', 'G8', 'G8',
  2005. 'G9', 'G9', 'G9', 'G9',
  2006. 'G28', 'G28', 'G28',
  2007. 'G29', 'G29', 'G29',
  2008. 'G37', 'G37', 'G37', 'G37'
  2009. ]
  2011. stress_animal = ['G7', 'G8', 'G9', 'G28', 'G29', 'G37', 'G37']
  2013. def build_df_two_clusters(animal_list, y1, y2, condition):
  2014. y1 = np.asarray(y1, dtype=float)
  2015. y2 = np.asarray(y2, dtype=float)
  2016. animal_list = np.asarray(animal_list)
  2018. if len(animal_list) != len(y1) or len(animal_list) != len(y2):
  2019. raise ValueError(
  2020. f"Length mismatch: animals={len(animal_list)}, y1={len(y1)}, y2={len(y2)}"
  2021. )
  2023. df = pd.DataFrame({
  2024. "Animal_id": animal_list,
  2025. "y1": y1, # cluster 1 size
  2026. "y2": y2 # cluster 2 size
  2027. })
  2029. df["condition"] = condition
  2030. df["repeat"] = df.groupby("Animal_id").cumcount() + 1
  2031. return df[["Animal_id", "condition", "repeat", "y1", "y2"]]
  2034. # --- Build condition-specific dataframes ---
  2035. df_stim = build_df_two_clusters(stim_animal, flat_cluster1stim, flat_cluster2stim, "stim")
  2036. df_stress = build_df_two_clusters(stress_animal, flat_cluster1stress, flat_cluster2stress, "stress")
  2038. # --- Merge into one dataframe ---
  2039. df = pd.concat([df_stim, df_stress], ignore_index=True)
  2040. df = df.sort_values(["condition", "Animal_id", "repeat"]).reset_index(drop=True)
  2042. print(df)
  2044. # %%
  2045. # df has columns: Animal_id, condition, repeat, y1, y2
  2046. df_long = (
  2047. df[["Animal_id", "condition", "repeat", "y1", "y2"]]
  2048. .melt(
  2049. id_vars=["Animal_id", "condition", "repeat"],
  2050. value_vars=["y1", "y2"],
  2051. var_name="cluster",
  2052. value_name="y"
  2053. )
  2054. )
  2056. # clean cluster labels
  2057. df_long["cluster"] = df_long["cluster"].str.strip().map({"y1": "cluster1", "y2": "cluster2"})
  2059. # drop any rows where y is missing (should usually be none)
  2060. df_long = df_long.dropna(subset=["y"]).reset_index(drop=True)
  2062. print(df_long.head(10))
  2063. print(df_long.columns)
  2065. # %%
  2066. # Fit mixed-effects model with interaction between cluster and condition
  2067. m_nested = smf.mixedlm(
  2068. "y ~ cluster * condition",
  2069. df_long,
  2070. groups=df_long["Animal_id"],
  2071. ).fit(reml=False, method="nm", maxiter=2000, disp=True)
  2073. print(m_nested.summary())
  2075. # %%
  2076. # Compute marginal R²: variance explained by fixed effects only
  2077. var_fixed = np.var(m_nested.fittedvalues)
  2078. var_random = m_nested.cov_re.iloc[0, 0]
  2079. var_resid = m_nested.scale
  2081. R2_marginal = var_fixed / (var_fixed + var_random + var_resid)
  2082. print(R2_marginal)
  2083. # Extract p value for the fixed effect of cluster identity (cluster2 vs reference)
  2084. pvals_all_clusters = m_nested.pvalues.get("cluster[T.cluster2]")
  2085. print(pvals_all_clusters)
  2089. # %%
  2090. # Subset to stimulation only
  2091. df_stim_only = df_long[df_long["condition"] == "stim"].copy()
  2092. # Fit mixed-effects model: cluster effect within stim
  2093. m_stim = smf.mixedlm(
  2094. "y ~ cluster",
  2095. df_stim_only,
  2096. groups=df_stim_only["Animal_id"],
  2097. re_formula="~1",
  2098. ).fit(reml=False, method="nm", maxiter=2000, disp=True)
  2100. print(m_stim.summary())
  2102. # %%
  2103. # Compute marginal R² for the stimulation model (fixed effects only)
  2104. var_fixed = np.var(m_stim.fittedvalues)
  2105. var_random = m_stim.cov_re.iloc[0, 0]
  2106. var_resid = m_stim.scale
  2108. R2_marginal = var_fixed / (var_fixed + var_random + var_resid)
  2109. print(R2_marginal)
  2111. # Extract p value for the cluster effect under stimulation
  2112. pvals_stim_clusters = m_stim.pvalues.get("cluster[T.cluster2]")
  2113. print(pvals_stim_clusters)
  2116. # %%
  2117. # Subset to stress only
  2118. df_stress_only = df_long[df_long["condition"] == "stress"].copy()
  2119. # Fit mixed-effects model: cluster effect within stress
  2120. m_stress = smf.mixedlm(
  2121. "y ~ cluster",
  2122. df_stress_only,
  2123. groups=df_stress_only["Animal_id"],
  2124. re_formula="~1",
  2125. ).fit(reml=False, method="nm", maxiter=2000, disp=True)
  2127. print(m_stress.summary())
  2129. # %%
  2130. # Compute marginal R² for the stress model (fixed effects only)
  2131. var_fixed = np.var(m_stress.fittedvalues)
  2132. var_random = m_stress.cov_re.iloc[0, 0]
  2133. var_resid = m_stress.scale
  2135. R2_marginal = var_fixed / (var_fixed + var_random + var_resid)
  2136. print(R2_marginal)
  2138. # Extract p value for the cluster effect under stress
  2139. pvals_stress_clusters = m_stress.pvalues.get("cluster[T.cluster2]")
  2140. print(pvals_stress_clusters)
  2143. # %%
  2144. # Extract cluster category labels from column headers
  2145. categories_percent_hier = cell_percentage_hier.columns.tolist()
  2147. # Extract total number of GABAergic cells across trials
  2148. cluster1_cells_allanimals_stim = cell_percentage_hier['Stimulation (#1)']
  2149. cluster2_cells_allanimals_stim = cell_percentage_hier['Stimulation (#2)']
  2151. cluster1_cells_allanimals_stress = cell_percentage_hier['Stress (#1)']
  2152. cluster2_cells_allanimals_stress = cell_percentage_hier['Stress (#2)']
  2154. # Remove all the NaNs
  2155. cluster1_cells_allanimals_stim = cluster1_cells_allanimals_stim[~np.isnan(cluster1_cells_allanimals_stim)]
  2156. cluster2_cells_allanimals_stim = cluster2_cells_allanimals_stim[~np.isnan(cluster2_cells_allanimals_stim)]
  2158. cluster1_cells_allanimals_stress = cluster1_cells_allanimals_stress[~np.isnan(cluster1_cells_allanimals_stress)]
  2159. cluster2_cells_allanimals_stress = cluster2_cells_allanimals_stress[~np.isnan(cluster2_cells_allanimals_stress)]
  2162. # %%
  2163. # Set up figure layout
  2164. fig = plt.figure(figsize=(10, 15)) # Define overall figure size
  2165. gs = gridspec.GridSpec(8, 3, figure=fig) # Create an 8x3 grid for subplot placement
  2167. # Sort cell indices by hierarchical clustering labels (stimulus condition)
  2168. idx = np.argsort(labels_hier_stim_2)
  2170. # 🔳 Plot 1: Heatmap of pairwise correlations (stimulus condition)
  2171. ax = fig.add_subplot(gs[0:2, 0]) # Top-left subplot
  2172. ax = sns.heatmap(max_corr_stim_2[idx, :][:, idx], # Correlation matrix sorted by cluster
  2173. cbar_kws={'label': 'Corr', 'location': 'top'}, # Colourbar settings
  2174. cmap="viridis", vmin=-1, vmax=1) # Colourmap and value range
  2175. ax.set(xlabel='cell ID', ylabel='cell ID') # Axis labels
  2176. ax.set_xticks([_ for _ in range(n_cells_stim) if _ % 10 == 0]) # X-axis ticks every 10 cells
  2177. ax.set_xticklabels([_ for _ in range(n_cells_stim) if _ % 10 == 0]) # X-axis tick labels
  2178. ax.set_yticks([_ for _ in range(n_cells_stim) if _ % 10 == 0]) # Y-axis ticks every 10 cells
  2179. ax.set_yticklabels([_ for _ in range(n_cells_stim) if _ % 10 == 0]) # Y-axis tick labels
  2180. ax.invert_yaxis() # Show first cell at top
  2182. # Add plot outline
  2183. for spine in ax.spines.values():
  2184. spine.set_visible(True)
  2185. spine.set_linewidth(1.5)
  2186. spine.set_edgecolor('gray')
  2188. # Add colourbar outline
  2189. cbar = ax.collections[0].colorbar
  2190. cbar.outline.set_visible(True)
  2191. cbar.outline.set_edgecolor('grey')
  2192. cbar.outline.set_linewidth(1.0)
  2194. # 🔳 Plot 2: Heatmap of normalised fluorescence traces (stimulus condition)
  2195. ax = plt.subplot(gs[0:2, 1:3]) # Top-right subplot
  2196. im = ax.imshow(minmax_scale(scale(data_stim_2[idx, :], axis=1), axis=1), aspect='auto',
  2197. cmap="viridis", interpolation='none', vmin=0, vmax=1, origin='upper') # Heatmap of scaled traces
  2198. ax.set_xlabel('')
  2199. ax.set_xticks([])
  2200. ax.set_yticks(range(0, 40, 30))
  2201. ax.set_yticklabels([])
  2202. ax.invert_yaxis()
  2204. # Add colourbar for fluorescence heatmap
  2205. cbar = plt.colorbar(im, ax=ax, orientation='horizontal', location='top', aspect=50, pad=0.05)
  2206. cbar.set_label('$\Delta$ F$/$F$_0$ (scaled)')
  2208. # 🔳 Plot 3: Average traces per cluster (stimulus condition)
  2209. ax = plt.subplot(gs[2, 1:3]) # Middle-right subplot
  2210. for Ki in range(Ks_2): # Loop through clusters
  2211. js = np.where(labels_hier_stim_2 == Ki) # Get cell indices for cluster Ki
  2212. avg = np.mean(np.squeeze(data_stim_2[js, :]), axis=0) # Average trace
  2213. avg = minmax_scale(avg) # Normalise trace
  2214. ax.plot(avg, c=GABA_greens[Ki-1]) # Plot with cluster colour
  2215. ax.set_xlim(0, 120)
  2216. ax.set_xlabel('time [sec]')
  2217. ax.set_xticks(range(0, 1400, 200), labels=range(0, 140, 20))
  2218. ax.set_ylabel('$\Delta F/F_0$')
  2219. ax.set_yticks([0, 0.5, 1])
  2221. # Sort cell indices by hierarchical clustering labels (stress condition)
  2222. idx = np.argsort(labels_hier_stress_2)
  2224. # 🔳 Plot 4: Heatmap of pairwise correlations (stress condition)
  2225. ax = fig.add_subplot(gs[3:5, 0]) # Middle-left subplot
  2226. ax = sns.heatmap(max_corr_stress_2[idx, :][:, idx],
  2227. cbar_kws={'label': 'Corr', 'location': 'top'},
  2228. cmap="viridis", vmin=-1, vmax=1)
  2229. ax.set(xlabel='cell ID', ylabel='cell ID')
  2230. ax.set_xticks([_ for _ in range(n_cells_stress) if _ % 20 == 0])
  2231. ax.set_xticklabels([_ for _ in range(n_cells_stress) if _ % 20 == 0])
  2232. ax.set_yticks([_ for _ in range(n_cells_stress) if _ % 20 == 0])
  2233. ax.set_yticklabels([_ for _ in range(n_cells_stress) if _ % 20 == 0])
  2234. ax.invert_yaxis()
  2236. # Add plot outline
  2237. for spine in ax.spines.values():
  2238. spine.set_visible(True)
  2239. spine.set_linewidth(1.5)
  2240. spine.set_edgecolor('gray')
  2242. # Add colourbar outline
  2243. cbar = ax.collections[0].colorbar
  2244. cbar.outline.set_visible(True)
  2245. cbar.outline.set_edgecolor('grey')
  2246. cbar.outline.set_linewidth(1.0)
  2248. # 🔳 Plot 5: Heatmap of normalised fluorescence traces (stress condition)
  2249. ax = plt.subplot(gs[3:5, 1:3]) # Middle-right subplot
  2250. im = ax.imshow(minmax_scale(scale(data_stress_2[idx, :], axis=1), axis=1), aspect='auto',
  2251. cmap="viridis", interpolation='none', vmin=0, vmax=1, origin='upper')
  2252. ax.set_xlabel('')
  2253. ax.set_xticks([])
  2254. ax.set_yticks(range(0, 40, 30))
  2255. ax.set_yticklabels([])
  2256. ax.invert_yaxis()
  2258. # Add colourbar for fluorescence heatmap
  2259. cbar = plt.colorbar(im, ax=ax, orientation='horizontal', location='top', aspect=50, pad=0.05)
  2260. cbar.set_label('$\Delta$ F$/$F$_0$ (scaled)')
  2262. # 🔳 Plot 6: Average traces per cluster (stress condition)
  2263. ax = plt.subplot(gs[5, 1:3]) # Bottom-right subplot
  2264. for Ki in range(Ks_2):
  2265. js = np.where(labels_hier_stress_2 == Ki)
  2266. avg = np.mean(np.squeeze(data_stress_2[js, :]), axis=0)
  2267. avg = minmax_scale(avg)
  2268. ax.plot(avg, c=GABA_greens[Ki-1])
  2269. ax.set_xlim(0, 120)
  2270. ax.set_xlabel('time [sec]')
  2271. ax.set_xticks(range(0, 1500, 300), labels=range(0, 150, 30))
  2272. ax.set_ylabel('$\Delta F/F_0$')
  2273. ax.set_yticks([0, 0.5, 1])
  2275. # 🔳 Plot 7: Violin plot of cell percentages (stimulation condition, cluster group A)
  2276. ax = plt.subplot(gs[6:8, 0]) # Bottom-left subplot
  2277. perc1 = cluster1_cells_allanimals_stim
  2278. df = pd.DataFrame(dict(x=np.repeat([0], len(perc1)), y=perc1))
  2279. sns.violinplot(x="x", y="y", data=df, order=np.arange(2), inner='box', linecolor='k', color=eff_green)
  2281. # perc2 = cell_percentage_hier[categories_percent_hier[8]] / cell_percentage_hier[categories_percent_hier[1]] * 100
  2282. perc2 =cluster2_cells_allanimals_stim
  2283. df = pd.DataFrame(dict(x=np.repeat([1], len(perc2)), y=perc2))
  2284. sns.violinplot(x="x", y="y", data=df, order=np.arange(2), inner='box', linecolor='k', color=inter_green)
  2286. ax.set_title(categories_percent_hier[7][:-5])
  2287. ax.set_ylim(-20, 140)
  2288. ax.set_xlabel('')
  2289. ax.set_ylabel('% of cells')
  2290. ax.set_xticks(range(2), labels=['cluster 1', 'cluster 2'])
  2291. ax.set_yticks(range(0, 150, 50), labels=range(0, 150, 50))
  2293. # Statistical tests between cluster groups
  2294. ax.plot([0, 1], [123, 123], 'k-', lw=1)
  2295. ax.plot([0, 0], [70, 123], 'k-', lw=1)
  2296. ax.plot([1, 1], [118, 123], 'k-', lw=1)
  2297. MW_p_value, _ = mannwhitneyu(perc1.dropna(), perc2.dropna())
  2298. KS_p_value = pvals_stim_clusters
  2299. ax.annotate(pval_to_star(KS_p_value), (0.5, 125), ha='center')
  2301. # 🔳 Plot 8: Violin plot of cell percentages (stress condition, cluster group B)
  2302. ax = plt.subplot(gs[6:8, 1]) # Bottom-center subplot
  2303. perc1 = cluster1_cells_allanimals_stress
  2304. df = pd.DataFrame(dict(x=np.repeat([0], len(perc1)), y=perc1))
  2305. sns.violinplot(x="x", y="y", data=df, order=np.arange(2), inner='box', linecolor='k', color=eff_green)
  2307. # Extract and normalise cell percentages for cluster 2
  2308. perc2 = cluster2_cells_allanimals_stress
  2309. df = pd.DataFrame(dict(x=np.repeat([1], len(perc2)), y=perc2))
  2310. sns.violinplot(x="x", y="y", data=df, order=np.arange(2), inner='box', linecolor='k', color=inter_green)
  2312. # Configure plot appearance
  2313. ax.set_title(categories_percent_hier[11][:-5]) # Remove suffix from category name
  2314. ax.set_ylim(-20, 140) # Set y-axis limits
  2315. ax.set_xlabel('')
  2316. ax.set_ylabel('% of cells')
  2317. ax.set_xticks(range(2), labels=['cluster 1', 'cluster 2'])
  2318. ax.set_yticks(range(0, 150, 50), labels=range(0, 150, 50))
  2320. # Statistical tests between cluster groups
  2321. ax.plot([0, 1], [123, 123], 'k-', lw=1) # Horizontal line
  2322. ax.plot([0, 0], [70, 123], 'k-', lw=1) # Left vertical line
  2323. ax.plot([1, 1], [118, 123], 'k-', lw=1) # Right vertical line
  2324. MW_p_value, _ = mannwhitneyu(perc1.dropna(), perc2.dropna()) # Mann-Whitney U test
  2325. KS_p_value = pvals_stress_clusters # Kolmogorov-Smirnov test
  2326. ax.annotate(pval_to_star(KS_p_value), (0.5, 125), ha='center') # Annotate significance
  2328. # Optional: Add final layout adjustments and save figure
  2329. fig.tight_layout()
  2330. # plt.savefig('cluster_summary_figure.svg', dpi=700) # Uncomment to save
  2331. plt.show()
  2333. # %% [markdown]
  2334. # # 🖼️ Figure S3 (cluster size)
  2336. # %%
  2337. # Create figure for cluster percentage comparison across all animals
  2338. fig = plt.figure(figsize=(8, 5))
  2340. ax = plt.subplot(1, 1, 1)
  2341. ax.tick_params(axis='both', length=0)
  2343. # Extract cell percentages for the two clusters
  2344. perc1 = cell_percentage_hier['All Animals (#1)']
  2345. perc2 = cell_percentage_hier['All Animals (#2)']
  2347. # Violin plot for cluster 1
  2348. df1 = pd.DataFrame(dict(x=np.repeat([0], len(perc1)), y=perc1))
  2349. sns.violinplot(x="x", y="y", data=df1, order=np.arange(2), color=eff_green)
  2351. # Violin plot for cluster 2
  2352. df2 = pd.DataFrame(dict(x=np.repeat([1], len(perc2)), y=perc2))
  2353. sns.violinplot(x="x", y="y", data=df2, order=np.arange(2), color=inter_green)
  2355. # Axis formatting and labels
  2356. ax.set_title(categories_percent_hier[4*0 + 3][:-5])
  2357. ax.set_ylim(-25, 150)
  2358. ax.set_xlabel('')
  2359. ax.set_ylabel('% of cells')
  2360. ax.set_xticks([0, 1], labels=['cluster 1', 'cluster 2'])
  2362. # Annotate p value for cluster comparison
  2363. ax.annotate(f'p={pvals_all_clusters:0.3e}', (0.5, 125), ha='center')
  2365. fig.tight_layout()
  2366. # plt.savefig('cluster_summary_figure_supplement.svg', dpi=700)
  2367. plt.show()
  2370. # %% [markdown]
  2371. # # 🖼️ Figure S4 (corrcoef)
  2373. # %%
  2374. # List of all relevant CSV file paths for UCN-opsin & GABA-GCaMP experiments
  2375. all_file_paths = [
  2376. './data/calcium_imaging/UCN3_stimulation/G7 20231110 10hz 222sti.csv', #1
  2377. './data/calcium_imaging/UCN3_stimulation/G7 20240201 10hz 223sti.csv', #2
  2378. './data/calcium_imaging/UCN3_stimulation/G8 20240220 10hz 223sti.csv', #3
  2379. './data/calcium_imaging/UCN3_stimulation/G8 20240301 10hz 223sti.csv', #4
  2380. './data/calcium_imaging/UCN3_stimulation/G8 20240202 10hz 223sti.csv', #5
  2381. './data/calcium_imaging/UCN3_stimulation/G8 20231218 10hz 222sti.csv', #6
  2382. './data/calcium_imaging/UCN3_stimulation/G9 20231218 10hz 222sti.csv', #7
  2383. './data/calcium_imaging/UCN3_stimulation/G9 20240202 10hz 223sti.csv', #8
  2384. './data/calcium_imaging/UCN3_stimulation/G9 20240307 10hz 223sti.csv', #9
  2385. './data/calcium_imaging/UCN3_stimulation/G9 20240228 10hz 223sti.csv', #10
  2386. './data/calcium_imaging/UCN3_stimulation/G28 20240709 10hz 223sti.csv', #11
  2387. './data/calcium_imaging/UCN3_stimulation/G28 20240701 10hz 223sti.csv', #12
  2388. './data/calcium_imaging/UCN3_stimulation/G28 20240704 10hz 223sti.csv', #13
  2389. './data/calcium_imaging/UCN3_stimulation/G29 20240703 10hz 223sti.csv', #14
  2390. './data/calcium_imaging/UCN3_stimulation/G29 20240711 10hz 223sti.csv', #15
  2391. './data/calcium_imaging/UCN3_stimulation/G29 20240705 10hz 223sti.csv', #16
  2392. './data/calcium_imaging/UCN3_stimulation/G37 10hz 223sti 20240720 A.csv', #17
  2393. './data/calcium_imaging/UCN3_stimulation/G37 10hz 223sti 20240720 B.csv', #18
  2394. './data/calcium_imaging/UCN3_stimulation/G37 10hz 223sti 20240721 A.csv', #19
  2395. './data/calcium_imaging/UCN3_stimulation/G37 10hz 223sti 20240721 B.csv', #20
  2396. './data/calcium_imaging/UCN3_stimulation/G7 20231218 334restraint.csv', #21
  2397. './data/calcium_imaging/UCN3_stimulation/G8 20231218 343restraint.csv', #22
  2398. './data/calcium_imaging/UCN3_stimulation/G9 20231218 334restraint.csv', #23
  2399. './data/calcium_imaging/restrain_stress/G28 333restraint 20240630.csv', #24
  2400. './data/calcium_imaging/restrain_stress/G29 333restraint 20240701.csv', #25
  2401. './data/calcium_imaging/restrain_stress/G37 333restraint 20240719.csv', #26
  2402. './data/calcium_imaging/restrain_stress/G37 333restraint 20240725.csv', #27
  2403. ]
  2406. # %%
  2407. # Selected stimulation files for bootstrapping analysis
  2408. # Each entry corresponds to a specific animal and session
  2410. boot_file_paths = [
  2411. all_file_paths[1], # G7 — 2024-02-01 — 223sti
  2412. all_file_paths[5], # G8 — 2024-03-01 — 223sti
  2413. all_file_paths[8], # G9 — 2024-03-07 — 223sti
  2414. all_file_paths[12], # G28 — 2024-07-09 — 223sti
  2415. all_file_paths[15], # G29 — 2024-07-11 — 223sti
  2416. all_file_paths[17] # G37 — 2024-07-20 A — 223sti
  2417. ]
  2419. # %%
  2420. # Create a 3x3 grid layout for plotting silhouette score distributions
  2421. plt.figure(figsize=(15, 15)) # Set overall figure size
  2423. # Define subplot axes for each animal/group
  2424. ax1 = plt.subplot(3, 3, 1)
  2425. ax2 = plt.subplot(3, 3, 2)
  2426. ax3 = plt.subplot(3, 3, 3)
  2427. ax4 = plt.subplot(3, 3, 4)
  2428. ax5 = plt.subplot(3, 3, 5)
  2429. ax6 = plt.subplot(3, 3, 6)
  2431. # Iterate over bootstrapped stimulation datasets
  2432. for ind, stim_csv_fname in enumerate(boot_file_paths):
  2433. # Load stimulation data and raw traces
  2434. data_stim, cells_traces_tot_old = load_data(stim_csv_fname)
  2436. # Filter for UCN3+ GABAergic neurons
  2437. data_stim_tot, n_cells = process_data(data_stim, cells_traces_tot_old)
  2439. # Normalise each trace to its own min/max (row-wise)
  2440. normalised_data_stim_tot = minmax_scale(data_stim_tot, axis=1)
  2442. # Slice stimulation period from normalised traces
  2443. normalised_data_stim_2 = normalised_data_stim_tot[:, stim_onset:stim_offset]
  2445. # Compute pairwise correlation-based distance matrix
  2446. # Using 1 - Pearson correlation as a dissimilarity metric
  2447. dist_cov_2 = pdist(scale(normalised_data_stim_2, axis=1),
  2448. lambda x, y: 1 - np.corrcoef(x, y)[0, 1])
  2450. # Initialise array to store silhouette scores across bootstraps
  2451. sil_score_2 = np.zeros((max_no_clusters - 1, n_resampling))
  2453. # Loop over cluster numbers and bootstrap iterations
  2454. for i, Ks_2 in enumerate(range(2, max_no_clusters + 1)):
  2455. for j in range(n_resampling):
  2456. # Resample 75% of cells for bootstrapping
  2457. sampled_ids = np.random.choice(range(n_cells), int(0.75 * n_cells), replace=False)
  2458. sampled_data = squareform(dist_cov_2)[sampled_ids, :][:, sampled_ids]
  2460. # Apply hierarchical clustering with average linkage
  2461. ac = AgglomerativeClustering(n_clusters=Ks_2, metric='precomputed',
  2462. linkage='average', compute_distances=True).fit(sampled_data)
  2464. # Compute silhouette score using precomputed distances
  2465. sil_score_2[i, j] = silhouette_score(sampled_data, ac.labels_, metric='precomputed')
  2467. # Select subplot axis based on file name prefix
  2468. file_name = os.path.split(stim_csv_fname)[-1]
  2469. ax = {
  2470. 'G7': ax1, 'G8': ax2, 'G9': ax3,
  2471. 'G28': ax4, 'G29': ax5, 'G37': ax6
  2472. }.get(file_name.split(' ')[0], None)
  2474. if ax:
  2475. # Plot silhouette score distributions as violin plots
  2476. for k, colour in enumerate(['red', 'green', 'blue']): # Cluster sizes: 2, 3, 4
  2477. df = pd.DataFrame({'x': np.repeat([k], len(sil_score_2[k, :])), 'y': sil_score_2[k, :]})
  2478. sns.violinplot(x="x", y="y", data=df, color=colour, linecolor='k', ax=ax)
  2480. # Perform Dunn's test for post-hoc comparisons
  2481. p_vals = posthoc_dunn(sil_score_2, p_adjust='fdr_bh')
  2483. # Annotate significance between cluster numbers
  2484. ax.plot([0, 1], [0.79, 0.79], 'k-') # 2 vs 3
  2485. ax.plot([1, 2], [0.7, 0.7], 'k-') # 3 vs 4
  2486. ax.plot([0, 2], [0.89, 0.89], 'k-') # 2 vs 4
  2487. ax.annotate(pval_to_star(p_vals[1][2]), (0.5, 0.8), ha='center')
  2488. ax.annotate(pval_to_star(p_vals[2][3]), (1.5, 0.71), ha='center')
  2489. ax.annotate(pval_to_star(p_vals[1][3]), (1, 0.9), ha='center')
  2491. # Configure subplot titles and axis labels
  2492. for idx, title in enumerate(['G7', 'G8', 'G9', 'G28', 'G29', 'G37'], start=1):
  2493. ax = plt.subplot(3, 3, idx)
  2494. ax.set_title(title) # Animal/group ID
  2495. ax.set_ylim(0, 1) # Score range
  2496. ax.set_ylabel('Silhouette score')
  2497. ax.set_xlabel('Number of Clusters')
  2498. ax.set_xticklabels([2, 3, 4, '', 2, 3, 4]) # Cluster labels
  2500. # Final layout and export
  2501. plt.tight_layout() # Optimise spacing
  2502. # plt.savefig('./gdrive/My Drive/FigS3_corr_revised.svg', dpi=300) # Save figure
  2503. plt.show() # Display plot
  2505. # %% [markdown]
  2506. # # 🖼️ Figure S5 (SLxCorr)
  2508. # %%
  2509. # Create a 3x3 grid layout for plotting silhouette score distributions
  2510. plt.figure(figsize=(15, 15)) # Set overall figure size
  2512. # Define subplot axes for each animal/group
  2513. ax1 = plt.subplot(3, 3, 1)
  2514. ax2 = plt.subplot(3, 3, 2)
  2515. ax3 = plt.subplot(3, 3, 3)
  2516. ax4 = plt.subplot(3, 3, 4)
  2517. ax5 = plt.subplot(3, 3, 5)
  2518. ax6 = plt.subplot(3, 3, 6)
  2520. # Iterate over bootstrapped stimulation datasets
  2521. for ind, stim_csv_fname in enumerate(boot_file_paths):
  2522. # Load stimulation data and raw traces
  2523. data_stim, cells_traces_tot_old = load_data(stim_csv_fname)
  2525. # Filter for UCN3+ GABAergic neurons
  2526. data_stim_tot, n_cells = process_data(data_stim, cells_traces_tot_old)
  2528. # Normalise each trace to its own min/max (row-wise)
  2529. normalised_data_stim_tot = minmax_scale(data_stim_tot, axis=1)
  2531. # Slice stimulation period from normalised traces
  2532. normalised_data_stim_2 = normalised_data_stim_tot[:, stim_onset:stim_offset]
  2534. # Compute pairwise signed lagged correlation distance matrix
  2535. dist_corr_2 = pdist(scale(normalised_data_stim_2, axis=1),
  2536. lambda x, y: 1 - signed_lagged_corr(x, y)[0])
  2538. # Initialise array to store silhouette scores across bootstraps
  2539. sil_score_2 = np.zeros((max_no_clusters - 1, n_resampling))
  2541. # Loop over cluster numbers and bootstrap iterations
  2542. for i, Ks_2 in enumerate(range(2, max_no_clusters + 1)):
  2543. for j in range(n_resampling):
  2544. # Resample 75% of cells for bootstrapping
  2545. sampled_ids = np.random.choice(range(n_cells), int(0.75 * n_cells), replace=False)
  2546. sampled_data = squareform(dist_corr_2)[sampled_ids, :][:, sampled_ids]
  2548. # Apply hierarchical clustering with average linkage
  2549. ac = AgglomerativeClustering(n_clusters=Ks_2, metric='precomputed',
  2550. linkage='average', compute_distances=True).fit(sampled_data)
  2552. # Compute silhouette score using precomputed distances
  2553. sil_score_2[i, j] = silhouette_score(sampled_data, ac.labels_, metric='precomputed')
  2555. # Select subplot axis based on file name prefix
  2556. file_name = os.path.split(stim_csv_fname)[-1]
  2557. ax = {
  2558. 'G7': ax1, 'G8': ax2, 'G9': ax3,
  2559. 'G28': ax4, 'G29': ax5, 'G37': ax6
  2560. }.get(file_name.split(' ')[0], None)
  2562. if ax:
  2563. # Plot silhouette score distributions as violin plots
  2564. for k, colour in enumerate(['red', 'green', 'blue']): # Cluster sizes: 2, 3, 4
  2565. df = pd.DataFrame({'x': np.repeat([k], len(sil_score_2[k, :])), 'y': sil_score_2[k, :]})
  2566. sns.violinplot(x="x", y="y", data=df, color=colour, linecolor='k', ax=ax)
  2568. # Perform Dunn's test for post-hoc comparisons
  2569. p_vals = posthoc_dunn(sil_score_2, p_adjust='fdr_bh')
  2571. # Annotate significance between cluster numbers
  2572. ax.plot([0, 1], [0.79, 0.79], 'k-') # 2 vs 3
  2573. ax.plot([1, 2], [0.7, 0.7], 'k-') # 3 vs 4
  2574. ax.plot([0, 2], [0.89, 0.89], 'k-') # 2 vs 4
  2575. ax.annotate(pval_to_star(p_vals[1][2]), (0.5, 0.8), ha='center')
  2576. ax.annotate(pval_to_star(p_vals[2][3]), (1.5, 0.71), ha='center')
  2577. ax.annotate(pval_to_star(p_vals[1][3]), (1, 0.9), ha='center')
  2579. # Configure subplot titles and axis labels
  2580. for idx, title in enumerate(['G7', 'G8', 'G9', 'G28', 'G29', 'G37'], start=1):
  2581. ax = plt.subplot(3, 3, idx)
  2582. ax.set_title(title) # Animal/group ID
  2583. ax.set_ylim(0, 1) # Score range
  2584. ax.set_ylabel('Silhouette score')
  2585. ax.set_xlabel('Number of Clusters')
  2586. ax.set_xticklabels([2, 3, 4, '', 2, 3, 4]) # Cluster labels
  2588. # Final layout and export
  2589. plt.tight_layout() # Optimise spacing
  2590. plt.show() # Display plot
  2592. # %% [markdown]
  2593. # # 🖼️ Figure S6 (kmeans)
  2595. # %%
  2596. # Create a 3x3 grid for plotting silhouette score distributions
  2597. plt.figure(figsize=(15, 15)) # Set overall figure size
  2599. # Define subplot axes for each animal/group
  2600. ax1 = plt.subplot(3, 3, 1)
  2601. ax2 = plt.subplot(3, 3, 2)
  2602. ax3 = plt.subplot(3, 3, 3)
  2603. ax4 = plt.subplot(3, 3, 4)
  2604. ax5 = plt.subplot(3, 3, 5)
  2605. ax6 = plt.subplot(3, 3, 6)
  2607. # Iterate over bootstrapped stimulation datasets
  2608. for ind, stim_csv_fname in enumerate(boot_file_paths):
  2609. # Load stimulation data and raw traces
  2610. data_stim, cells_traces_tot_old = load_data(stim_csv_fname)
  2612. # Filter for UCN3+ GABAergic neurons
  2613. data_stim_tot, n_cells = process_data(data_stim, cells_traces_tot_old)
  2615. # Normalise each trace to its own min/max (row-wise)
  2616. normalised_data_stim_tot = minmax_scale(data_stim_tot, axis=1)
  2618. # Slice stimulation period from normalised traces
  2619. normalised_data_stim_2 = normalised_data_stim_tot[:, stim_onset:stim_offset]
  2621. # Compute pairwise correlation-based distance matrix
  2622. dist_corr_2 = pdist(scale(normalised_data_stim_2, axis=1),
  2623. lambda x, y: 1 - np.corrcoef(x, y)[0, 1])
  2625. # Initialise array to store silhouette scores across bootstraps
  2626. sil_score_2 = np.zeros((max_no_clusters - 1, n_resampling))
  2628. # Loop over cluster numbers and bootstrap iterations
  2629. for i, Ks_2 in enumerate(range(2, max_no_clusters + 1)):
  2630. for j in range(n_resampling):
  2631. # Resample 75% of cells for bootstrapping
  2632. sampled_ids = np.random.choice(range(n_cells), int(0.75 * n_cells), replace=False)
  2633. sampled_data = squareform(dist_corr_2)[sampled_ids, :][:, sampled_ids]
  2635. # Apply K-means clustering to resampled data
  2636. km = KMeans(n_clusters=Ks_2, n_init=10, max_iter=100).fit(sampled_data)
  2638. # Compute silhouette score using precomputed distances
  2639. sil_score_2[i, j] = silhouette_score(sampled_data, km.labels_, metric='precomputed')
  2641. # Select subplot axis based on file name prefix
  2642. file_name = os.path.split(stim_csv_fname)[-1]
  2643. ax = {
  2644. 'G7': ax1, 'G8': ax2, 'G9': ax3,
  2645. 'G28': ax4, 'G29': ax5, 'G37': ax6
  2646. }.get(file_name.split(' ')[0], None)
  2648. if ax:
  2649. # Plot silhouette score distributions as violin plots
  2650. for k, colour in enumerate(['red', 'green', 'blue']):
  2651. df = pd.DataFrame({'x': np.repeat([k], len(sil_score_2[k, :])), 'y': sil_score_2[k, :]})
  2652. sns.violinplot(x="x", y="y", data=df, color=colour, linecolor='k', ax=ax)
  2654. # Perform Dunn's test for post-hoc comparisons
  2655. p_vals = posthoc_dunn(sil_score_2, p_adjust='fdr_bh')
  2657. # Annotate significance between cluster numbers
  2658. ax.plot([0, 1], [0.79, 0.79], 'k-') # 2 vs 3
  2659. ax.plot([1, 2], [0.7, 0.7], 'k-') # 3 vs 4
  2660. ax.plot([0, 2], [0.89, 0.89], 'k-') # 2 vs 4
  2661. ax.annotate(pval_to_star(p_vals[1][2]), (0.5, 0.8), ha='center')
  2662. ax.annotate(pval_to_star(p_vals[2][3]), (1.5, 0.71), ha='center')
  2663. ax.annotate(pval_to_star(p_vals[1][3]), (1, 0.9), ha='center')
  2665. # Configure subplot titles and axis labels
  2666. for idx, title in enumerate(['G7', 'G8', 'G9', 'G28', 'G29', 'G37'], start=1):
  2667. ax = plt.subplot(3, 3, idx)
  2668. ax.set_title(title) # Animal/group ID
  2669. ax.set_ylim(0, 1) # Score range
  2670. ax.set_ylabel('Silhouette score')
  2671. ax.set_xlabel('Number of Clusters')
  2672. ax.set_xticklabels([2, 3, 4, '', 2, 3, 4]) # Cluster labels
  2674. # Final layout and export
  2675. plt.tight_layout()
  2676. # plt.savefig('./gdrive/My Drive/FigS5_kmeans_revised.svg', dpi=300)
  2677. plt.show()
  2679. # %% [markdown]
  2680. # # 🖼️ Figure S7 (stim. frequency)
  2682. # %%
  2683. # Distance metrics for control condition at 5 Hz stimulation
  2684. # Each value likely represents a summary statistic (e.g., mean trace distance) per cell or ROI
  2685. control_dist_5Hz = [0.3880, 0.5961, 0.6274, 0.2502, 0.3371, 0.3458, 0.3078, 0.4944, 0.3302, 0.3775, 0.5336, 0.4705]
  2687. # Distance metrics for stimulated condition at 5 Hz
  2688. # Larger sample size suggests more cells or trials were analyzed under stimulation
  2689. stim_dist_5Hz = [0.3697, 0.3987, 0.3547, 0.4152, 0.3759, 0.3717, 0.3927, 0.4157, 0.3807, 0.3770,
  2690. 0.3625, 0.2734, 0.3256, 0.4289, 0.3702, 0.3674, 0.4469, 0.4589, 0.4747, 0.3439]
  2692. # Distance metrics for control condition at 20 Hz stimulation
  2693. # Values may reflect baseline or spontaneous activity patterns at higher frequency context
  2694. control_dist_20Hz = [0.6933, 0.5293, 0.5910, 0.4005, 0.3390, 0.4714, 0.5060, 0.4065, 0.4393, 0.5510, 0.4441, 0.4634]
  2696. # Distance metrics for stimulated condition at 20 Hz
  2697. # Used to assess how stimulation alters trace similarity or variability at higher frequency
  2698. stim_dist_20Hz = [0.3857, 0.4040, 0.3860, 0.3939, 0.4034, 0.3658, 0.4081, 0.4177, 0.4064, 0.3816,
  2699. 0.2990, 0.4429, 0.3006, 0.4532, 0.3107, 0.4206, 0.4544, 0.4467, 0.4376, 0.4959]
  2701. # %%
  2702. # =========================
  2703. # Helpers
  2704. # =========================
  2705. def make_palette(animals, colors):
  2706. uniq = pd.unique(animals)
  2707. if len(uniq) != len(colors):
  2708. raise ValueError("Number of animals and colors must match")
  2709. return dict(zip(uniq, colors))
  2711. palette_ctrl = make_palette(control_animal, colors_ctrl)
  2712. palette_stim = make_palette(stim_animal, mycolors)
  2714. x_order = [0, 1, 3, 4, 6, 7, 8, 9 ]
  2716. def violin_swarm(xpos, values, animals, palette, marker):
  2717. df = pd.DataFrame({
  2718. "x": np.repeat(xpos, len(values)),
  2719. "y": values,
  2720. "c": animals
  2721. })
  2723. sns.violinplot(
  2724. x="x", y="y", data=df,
  2725. order=x_order,
  2726. fill=False, color="k",
  2727. linewidth=1.5,
  2728. ax=ax
  2729. )
  2731. sns.swarmplot(
  2732. x="x", y="y", hue="c", data=df,
  2733. order=x_order,
  2734. palette=palette,
  2735. size=10,
  2736. marker=marker,
  2737. legend=False,
  2738. ax=ax
  2739. )
  2741. # =========================
  2742. # Figure
  2743. # =========================
  2744. fig = plt.figure(figsize=(10, 10))
  2745. ax = plt.subplot(1, 1, 1)
  2747. # --- 5 Hz ---
  2748. violin_swarm(0, control_dist_5Hz, control_animal, palette_ctrl, "s")
  2749. violin_swarm(1, stim_dist_5Hz, stim_animal, palette_stim, "o")
  2751. # --- 10 Hz ---
  2752. violin_swarm(4, control_dist, control_animal, palette_ctrl, "s")
  2753. violin_swarm(6, stim_dist, stim_animal, palette_stim, "o")
  2755. # --- 20 Hz ---
  2756. violin_swarm(8, control_dist_20Hz, control_animal, palette_ctrl, "s")
  2757. violin_swarm(9, stim_dist_20Hz, stim_animal, palette_stim, "o")
  2759. # =========================
  2760. # Axes
  2761. # =========================
  2762. ax.set_ylim(0.1, 1.2)
  2763. ax.set_ylabel("Normalised Riemannian Distance")
  2764. ax.set_xlabel("")
  2765. ax.set_xticks(
  2766. [0, 1, 3, 4, 6, 7],
  2767. labels=[
  2768. "5Hz\n Control", "5Hz\n Stim.",
  2769. "10Hz\n Control", "10Hz\n Stim.",
  2770. "20Hz\n Control", "20Hz\n Stim."
  2771. ]
  2772. )
  2774. # =========================
  2775. # Statistics
  2776. # =========================
  2777. # 5 Hz
  2778. ax.plot([0, 1], [0.85, 0.85], "k-", lw=1)
  2779. ax.plot([0, 0], [0.85, 0.84], "k-", lw=1)
  2780. ax.plot([1, 1], [0.85, 0.84], "k-", lw=1)
  2781. _, p = kruskal(control_dist_5Hz, stim_dist_5Hz)
  2782. ax.annotate(pval_to_star(p), (0.5, 0.86), ha="center")
  2784. # 10 Hz
  2785. ax.plot([3, 4], [1.1, 1.1], "k-", lw=1)
  2786. ax.plot([3, 3], [1.1, 1.09], "k-", lw=1)
  2787. ax.plot([4, 4], [1.1, 1.09], "k-", lw=1)
  2788. _, p = kruskal(control_dist, stim_dist)
  2789. ax.annotate(pval_to_star(p), (3.5, 1.11), ha="center")
  2791. # 20 Hz
  2792. ax.plot([6, 7], [0.85, 0.85], "k-", lw=1)
  2793. ax.plot([6, 6], [0.85, 0.84], "k-", lw=1)
  2794. ax.plot([7, 7], [0.85, 0.84], "k-", lw=1)
  2795. _, p = kruskal(control_dist_20Hz, stim_dist_20Hz)
  2796. ax.annotate(pval_to_star(p), (6.5, 0.86), ha="center")
  2798. plt.tight_layout()
  2799. # plt.savefig("Fig_dist_violin_swarm.svg", dpi=600, bbox_inches="tight")
  2800. plt.show()
  2802. # %% [markdown]
  2803. # # 📚 Used Libraries
  2805. # %%
  2806. # List all currently imported libraries in the global namespace
  2807. import types
  2809. def list_imported_modules():
  2810. return sorted(
  2811. val.__name__ for name, val in globals().items()
  2812. if isinstance(val, types.ModuleType)
  2813. )
  2815. # Display the list
  2816. set(list_imported_modules())

analysis_pipeline.ipynb at commit 538ba57, under CC0-1.0 · at the source

Overview

Authors: Junru Yu1,2, Saeed Farjami3,4,5, Kateryna Nechyporenko3,4, Xiao Feng Li1, Hafsa Yaseen1,6, Yanyan Lin1, Jinbin Ye1, Owen Hollings1, Ross de Burgh1, Baban Singh1, Kevin T O’Byrne1, Krasimira Tsaneva-Atanasova3,4,7, Margaritis Voliotis3,4
ORCID iDs: Junru Yu, Saeed Farjami, Kateryna Nechyporenko, Xiao Feng Li, Hafsa Yaseen, Yanyan Lin, Kevin T O’Byrne, Krasimira Tsaneva-Atanasova, Margaritis Voliotis
  1. Department of Women and Children’s Health, School of Life Course and Population Sciences, King’s College London, Guy’s Campus, London, UK
  2. Department of Rehabilitation Medicine, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, Zhejiang China
  3. Department of Mathematics and Statistics, University of Exeter, Exeter, UK
  4. Living Systems Institute, University of Exeter, Exeter, UK
  5. The Pirbright Institute, Pirbright, Surrey, UK
  6. Biological Sciences, University of Missouri, Columbia, MO USA
  7. EPSRC Hub for Quantitative Modelling in Healthcare, University of Exeter, Exeter, UK
Institutions: King's College London (United Kingdom); First Affiliated Hospital of Wenzhou Medical University (China); The Pirbright Institute (United Kingdom); University of Exeter (United Kingdom); University of Missouri (United States)
Journal: Nature communications, volume 17, issue 1, article 5690
Dates: received 9 February 2025; accepted 25 February 2026; published online 10 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-70364-9 · PMID 41807404 · PMCID PMC13319269 · OpenAlex W7134887411
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), systems (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, fMRI & imaging, Single-unit activity, calcium imaging
Keywords: Endocrine reproductive disorders, Neural circuits, Network models
MeSH: Amygdala*, GABAergic Neurons*, Reproduction*, Stress, Physiological*, Stress, Psychological*, Animals, Female, gamma-Aminobutyric Acid, Gonadotropin-Releasing Hormone, Hypothalamic-Pituitary-Gonadal Axis, Luteinizing Hormone, Mice, Mice, Inbred C57BL, Optogenetics, Urocortins (* major topic)
Topic: Neuroendocrine regulation and behavior (Social Psychology, Psychology), according to OpenAlex
Funding: Biotechnology and Biological Sciences Research Council (BB/W005883/1, BB/S019979/1, BB/W005913/1); RCUK | Engineering and Physical Sciences Research Council (EPSRC) (EP/T017856/)
Citations: cited by 1 paper (Europe PMC); 45 references in the paper
Research resources: i) coating antibody RRID:AB_2665514, ii) anti-LH antibody RRID:AB_2665533, RRID:AB_772206

Abstract

Stress can disrupt menstrual cycles, impair fertility and cause reproductive disfunction. The posterodorsal medial amygdala (MePD) integrates stress signals and regulates the gonadotropin-releasing hormone (GnRH) pulse generator through a dense network of GABA and Urocortin-3 (UCN3) neurons, yet the mechanisms underlying the circuitry remain poorly understood. Here, we combine in vivo mini-endoscopic calcium imaging, optogenetics, clustering analysis, and computational modeling to investigate the MePD circuitry in female mice. We uncover two anti-correlated GABA subpopulations in the MePD that are involved in the response to restraint stress and UCN3 neuron stimulation. Computational modeling suggests that mutual inhibition between these GABA groups drives their anti-correlated activity and predicts how these interactions shape downstream responses to stimulation of GABA and UCN3 neurons. In vivo optogenetics confirms that GABA neurons are critical for transmitting UCN3 signals to regulate luteinizing hormone (LH) pulse frequency. Together, our findings reveal amygdala GABAergic circuit mechanisms that mediate stress effects on reproductive health, linking emotional processing and neuroendocrine control.

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

Repository

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

mv-kr/MePD_GABA

License: CC0-1.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 538ba57379a0669c2cdbd2f5c34e6d5d9f37a1e9, 8 January 2026
Languages: MATLAB (15), Jupyter (7)
Size: 270 files, 22 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 7 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (7 files), pandas (6 files), SciPy (5 files), Matplotlib (4 files), statsmodels (4 files), Signal Processing Toolbox (2 files), seaborn (2 files), Pingouin (1 file), pyRiemann (1 file), scikit-learn (1 file), scikit-posthocs (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
24 files
At the source: github.com/mv-kr/MePD_GABA

Code availability

The code for reproducing the data analysis (Python™), mathematical modeling (MATLAB™), and figure generation for both the main text and Supplementary information, along with the datasets used in the analysis, is publicly available at Figshare45(10.24378/exe.31018072) or in Github repository (https://github.com/mv-kr/MePD_GABA).

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

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 22 scripts, each with its path and the digest of its content;
  • 17 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 main data supporting the results in this study are available within the paper and its Supplementary Information. Any additional requests for information can be directed to and will be fulfilled by the corresponding authors. Source data are provided with this paper. Source data are available at Figshare44 (10.24378/exe.31021447). Source data are provided with this paper for reproducing all Figures in the manuscript and Supplementary Information. Source data are provided with this paper.

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 1, 30 September 2026: the first record

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

Cite

This paper

Yu, J., Farjami, S., Nechyporenko, K., Li, X. F., Yaseen, H., Lin, Y., Ye, J., Hollings, O., de Burgh, R., Singh, B., O’Byrne, K. T., Tsaneva-Atanasova, K., & Voliotis, M. (2026). The role of amygdala GABA neurons in controlling stress and reproduction in female mice. Nature communications, 17(1), 5690. https://doi.org/10.1038/s41467-026-70364-9

BibTeX

@article{yu2026role,
author = {Yu, Junru and Farjami, Saeed and Nechyporenko, Kateryna and Li, Xiao Feng and Yaseen, Hafsa and Lin, Yanyan and Ye, Jinbin and Hollings, Owen and de Burgh, Ross and Singh, Baban and O’Byrne, Kevin T and Tsaneva-Atanasova, Krasimira and Voliotis, Margaritis},
title = {{The role of amygdala GABA neurons in controlling stress and reproduction in female mice}},
journal = {Nature communications},
year = {2026},
month = mar,
volume = {17},
number = {1},
pages = {5690},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-70364-9},
url = {https://doi.org/10.1038/s41467-026-70364-9},
pmid = {41807404},
pmcid = {PMC13319269}
}

RIS

TY - JOUR
AU - Yu, Junru
AU - Farjami, Saeed
AU - Nechyporenko, Kateryna
AU - Li, Xiao Feng
AU - Yaseen, Hafsa
AU - Lin, Yanyan
AU - Ye, Jinbin
AU - Hollings, Owen
AU - de Burgh, Ross
AU - Singh, Baban
AU - O’Byrne, Kevin T
AU - Tsaneva-Atanasova, Krasimira
AU - Voliotis, Margaritis
TI - The role of amygdala GABA neurons in controlling stress and reproduction in female mice
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/03/10
VL - 17
IS - 1
SP - 5690
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-70364-9
UR - https://doi.org/10.1038/s41467-026-70364-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-70364-9",
"type": "article-journal",
"title": "The role of amygdala GABA neurons in controlling stress and reproduction in female mice",
"container-title": "Nature communications",
"author": [
{
"family": "Yu",
"given": "Junru"
},
{
"family": "Farjami",
"given": "Saeed"
},
{
"family": "Nechyporenko",
"given": "Kateryna"
},
{
"family": "Li",
"given": "Xiao Feng"
},
{
"family": "Yaseen",
"given": "Hafsa"
},
{
"family": "Lin",
"given": "Yanyan"
},
{
"family": "Ye",
"given": "Jinbin"
},
{
"family": "Hollings",
"given": "Owen"
},
{
"family": "de Burgh",
"given": "Ross"
},
{
"family": "Singh",
"given": "Baban"
},
{
"family": "O’Byrne",
"given": "Kevin T"
},
{
"family": "Tsaneva-Atanasova",
"given": "Krasimira"
},
{
"family": "Voliotis",
"given": "Margaritis"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "5690",
"DOI": "10.1038/s41467-026-70364-9",
"PMID": "41807404",
"PMCID": "PMC13319269",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-70364-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
10
]
]
}
}

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.1016/j.isci.2026.117375 [code]
Motor priming is associated with widespread recruitment into neural ensembles and more rapid ensemble transitions.
Journal: iScience
In common: scikit-posthocs, Pingouin, Signal Processing Toolbox, 7 other tools, systems
[2] doi:10.7554/elife.100605 [code]
Age-related changes in ‘cortical’ 1/f dynamics are linked to cardiac activity
Journal: n/a
In common: pyRiemann, Pingouin, statsmodels, 6 other tools
[3] doi:10.1038/s41593-026-02362-5 [code]
Replay of procedural memory is independent of the hippocampus.
Journal: Nature neuroscience
In common: scikit-posthocs, Pingouin, statsmodels, 6 other tools, mouse
[4] doi:10.1038/s41467-026-74823-1 [code]
Cerebellar activity is triggered by reach endpoint during learning of a complex locomotor task.
Journal: Nature communications
In common: scikit-posthocs, Pingouin, statsmodels, 6 other tools, mouse
[5] doi:10.1038/s41467-026-75662-w [code]
Distinct Roles of Deep and Superficial Cortical Layers in Tone Prediction, Comparison, and Adaptation in Human Auditory Cortices.
Journal: Nature communications
In common: scikit-posthocs, Pingouin, statsmodels, 6 other tools
[6] doi:10.1186/s13321-026-01177-7 [code]
A pipeline for developing AI-driven models to predict molecular initiating events: a case study on neural tube defects.
Journal: Journal of cheminformatics
In common: scikit-posthocs, Pingouin, statsmodels, 6 other tools
[7] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: Pingouin, Signal Processing Toolbox, statsmodels, 6 other tools, systems, mouse
[8] doi:10.1002/advs.202518450 [code]
A Brain-Wide Atlas of Astrocytic Oxytocin Receptors Reveals a Glial Basis for Nucleus Accumbens Modulation of Affiliative Behavior.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: scikit-posthocs, Pingouin, statsmodels, 5 other tools, mouse
[9] doi:10.1038/s41586-026-10501-y [code]
Long-term editing of brain circuits using an engineered electrical synapse.
Journal: Nature
In common: Pingouin, Signal Processing Toolbox, statsmodels, 6 other tools, mouse
[10] doi:10.1371/journal.pbio.3003684 [code]
The retrieval of previously learned motor memories is facilitated by the reinstatement of default mode network manifold structures.
Journal: PLoS biology
In common: pyRiemann, Pingouin, seaborn, 5 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.