OSCR

Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.

Code ↔ Paper

11 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 11 matches
  1. [1] § Methods › Social maps ↔ umap/umap_.py, lines 1485–1737 · score 0.79 · low dimensional space, dimensional embedding, dimensionality reduction, UMAP embedded, uniformly, clustered
  2. [2] § Methods › Social maps ↔ social_behavior_umap.ipynb, lines 357–403 · score 0.78 · watershed segmentation, kernel density, UMAP embedded, stereotypic interaction, dimensional space, clustered
  3. [3] § Methods › Pose tracking and processing ↔ social_behavior_umap.ipynb, lines 177–241 · score 0.71 · relative orientation, relative angle, velocities, Lateral, speed, frames
  4. [4] § Results › Male- and female-directed courtship-like interactions exhibit different dynamics ↔ social_behavior_umap.ipynb, lines 357–403 · score 0.70 · kernel density, low density, stereotypical interaction, social modes, watershed, UMAP
  5. [5] § Methods › Maps of behavioral responses to song ↔ social_behavior_umap.ipynb, lines 177–241 · score 0.66 · wavelet transform, UMAP embedding, watershed, spectrogram, pose, signal
  6. [6] § Methods › Maps of behavioral responses to song ↔ umap/umap_.py, lines 1485–1737 · score 0.65 · dimensional space, UMAP embedding, manifold, uniformly, density, behavioral
  7. [7] § Methods › HMM-GLM modeling › Generalized linear model (GLM) ↔ demo/basics.ipynb, lines 100–127 · score 0.58 · raised cosine, small delays, broaden, GLM, transformation, filters
  8. [8] § Methods › Pose tracking and processing › Quantifying courter position during interactions ↔ social_behavior_umap.ipynb, lines 632–688 · score 0.57 · fly length, polar, histogram, angle, binned, courter
  9. [9] § Methods › Pose tracking and processing ↔ umap/spectral.py, lines 163–279 · score 0.55 · metrics function, relative positioning, Euclidean, vector, component, connecting
  10. [10] § Methods › HMM-GLM modeling ↔ ssm/transitions.py, lines 213–263 · score 0.53 · hidden Markov model, generalized linear model, transitioning, probability
  11. [11] § Results › Male- and female-directed courtship-like interactions exhibit different dynamics ↔ social_behavior_umap.ipynb, lines 841–898 · score 0.52 · male female, female directed, spent, edges, social, Ratio

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 · 985 lines · 41 KB · no license · 6 matches

  1. # %% [markdown]
  2. # # Social UMAP
  3. # This notebook contains the code to generate the social UMAPs from Nair et al. 2025, "Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior" (bioRxiv, May 2025).
  5. # %%
  6. # imports and configurations
  7. import os
  8. import pandas as pd
  9. import numpy as np
  10. from scipy import signal
  11. from sklearn.preprocessing import MinMaxScaler, StandardScaler
  12. from embedding_utils import (get_freq_scales, compute_wavelet,
  13. umap_embedding,
  14. time_delay_embedding,
  15. SpatialClustering,
  16. umap_multiple_fits,
  17. basis_transformation,
  18. kernel_density_estimation)
  19. from glm_utils import bases
  20. import matplotlib.pyplot as plt
  21. import seaborn as sns
  22. import matplotlib.colors as mcolors
  23. from scipy.stats import mannwhitneyu, binned_statistic_2d
  25. plt.style.use('ncb.mplstyle')
  27. colors = [mcolors.CSS4_COLORS["violet"], mcolors.CSS4_COLORS["royalblue"],
  28. mcolors.CSS4_COLORS["thistle"], mcolors.CSS4_COLORS["lightsteelblue"]]
  29. color_palette = sns.set_palette(sns.color_palette(colors))
  32. # %% [markdown]
  33. # ## How to arrange the data and specify parameters
  34. #
  35. # The data should be stored in the `dat` folder as one CSV file per trial. The columns of each CSV file are used as input features for the UMAP.
  36. #
  37. # The `trial_names` indicate the trials belonging to different groups, namely male-female and male-male interactions. The `features` list indicates the input features used to fit the UMAP, and `condition` indicates the condition on which the UMAPs are compared.
  38. #
  39. # - `start_frame_idx`: index of the first frame of each trial to be used in the analysis
  40. #
  41. # - `end_frame_idx`: index of the last frame of each trial to be used in the analysis
  42. #
  43. # - `confidence_threshold`: confidence level of tracked poses below which frames are ignored in the analysis
  44. #
  45. # ### Lowpass filtering
  46. # - `required_sampling_rate`: sampling rate of the data
  47. #
  48. # - `filter_features`: flag indicating whether to lowpass filter UMAP input features
  49. #
  50. # - `filter_cutoff`: cutoff frequency if lowpass filtering is used
  51. #
  52. # - `filter_ord`: filter order
  53. #
  54. # ### Normalization of input features
  55. # - `normalize_features`: whether to normalize the input features
  56. #
  57. # - `normalization`: what type of normalization is used, "standardize", "minmax", or `None`
  58. #
  59. # - `features_range`: for minmax normalization, the feature range to be used as minimum and maximum values. If used, this should be a dictionary whose keys are any or some of the feature names (`features`) and values are a list of two elements: the minimum and maximum value for the feature
  60. #
  61. # ### Time delay embedding parameters
  62. # - `time_delay_poses`: whether to use time-delay embedding of input features
  63. #
  64. # - `nb_delays`: if `time_delay_poses` is `True`, specify the number of delays to embed
  65. #
  66. # - `basis_transform`: whether to transform the basis vectors of time-delay embedded features
  67. #
  68. # - `basis_function`: if `basis_transform` is `True`, the function to be used for basis transformation
  69. #
  70. # ### Wavelet transform parameters
  71. # - `wavelet_transform`: whether to wavelet transform the input features. Either time-delay embedding or wavelet transform is used
  72. #
  73. # - `min_freq`: minimum frequency used for wavelet transform
  74. #
  75. # - `max_freq`: maximum frequency used for wavelet transform
  76. #
  77. # - `nfreqs`: number of frequency components to be used for wavelet transform
  78. #
  79. # - `freq_spacing`: frequency spacing to be used: 'dyadic', 'log', or 'uniform'
  80. #
  81. # ### UMAP parameters
  82. #
  83. # - `train_size`: ratio of data to be used for fitting the UMAP
  84. #
  85. # - `min_dist`: UMAP hyperparameter `min_dist` (https://umap-learn.readthedocs.io/en/latest/parameters.html)
  86. #
  87. # - `n_neighbors`: UMAP hyperparameter `n_neighbors` (https://umap-learn.readthedocs.io/en/latest/parameters.html)
  88. #
  89. # - `random_state`: random seed used to fit UMAP reproducibly
  90. #
  91. # - `nan_handling`: whether to 'remove' or 'interpolate' NaN values before analysis
  92. #
  93. # ### kernel density estimation (KDE) and watershed segmentation parameters
  94. # - `bw`: bandwidth to be used for KDE (https://kdepy.readthedocs.io/en/latest/introduction.html#Selecting-a-suitable-bandwidth)
  95. #
  96. # - `nb_gridpoints`: number of grid points on each 2D axis to be used for KDE
  97. #
  98. # - `watershed_threshold`: threshold to be used for watershed segmentation algorithm (see `watershed_segmentation` in `embedding_utils.py`)
  99. #
  100. # - `results_save_path`: path to save social UMAP results
  101. #
  102. # The notebook reads data from the CSV files and stores them in `preprocessed_data_path`.
  104. # %%
  105. data_save_path = "dat"
  107. preprocessed_data_path = os.path.join(data_save_path, "preprocessed_data.npz")
  109. # parameters to get data
  110. trial_names = {}
  111. trial_names['male-female'] = \
  112. ['localhost-20200706_122314',
  113. 'localhost-20200708_122011',
  114. 'localhost-20200707_113229',
  115. 'localhost-20200710_120433',
  116. 'localhost-20200710_134643',
  117. 'localhost-20200709_123809',
  118. 'localhost-20200709_115547',
  119. 'localhost-20200710_131755',
  120. 'localhost-20220214_114900',
  121. 'localhost-20220214_122614',
  122. 'localhost-20220214_130838',
  123. 'localhost-20220214_132737',
  124. 'localhost-20220216_093349',
  125. 'localhost-20220216_101942',
  126. 'localhost-20220218_114657',
  127. 'localhost-20220224_120331',
  128. 'localhost-20220225_121823',
  129. 'localhost-20220302_113245',
  130. 'localhost-20220302_125321',
  131. 'localhost-20220502_115227',
  132. 'localhost-20210127_145608',
  133. 'localhost-20210210_094512',
  134. 'localhost-20210210_100441',
  135. 'localhost-20210617_121041',
  136. 'localhost-20210617_123045',
  137. 'localhost-20220502_123531',
  138. 'localhost-20220506_130520',
  139. 'localhost-20220509_115048',
  140. 'localhost-20220509_124215',
  141. 'localhost-20220511_115515',
  142. ]
  144. trial_names['male-male'] = \
  145. ['localhost-20200706_132556',
  146. 'localhost-20200708_120051',
  147. 'localhost-20200812_143257',
  148. 'localhost-20200709_111319',
  149. 'localhost-20200812_155855',
  150. 'localhost-20210222_144954',
  151. 'localhost-20210222_152829',
  152. 'localhost-20210308_120623',
  153. 'localhost-20210308_122507',
  154. 'localhost-20220214_120801',
  155. 'localhost-20220214_135209',
  156. 'localhost-20220222_122109',
  157. 'localhost-20220223_121135',
  158. 'localhost-20220224_124055',
  159. 'localhost-20220224_130031',
  160. 'localhost-20220225_123838',
  161. 'localhost-20220225_125855',
  162. 'localhost-20220302_115211',
  163. 'localhost-20220302_131417',
  164. 'localhost-20220302_133433',
  165. 'localhost-20220502_125611',
  166. 'localhost-20210201_124220',
  167. 'localhost-20210201_132121',
  168. 'localhost-20210209_114709',
  169. 'localhost-20210209_122417',
  170. 'localhost-20220304_121729',
  171. 'localhost-20220304_123618',
  172. 'localhost-20220502_133927',
  173. 'localhost-20220509_121352',
  174. 'localhost-20220509_130538'
  175. ]
  177. features = [
  178. 'courter_velocity_forward',
  179. 'target_velocity_forward',
  180. 'courter_abs_velocity_lateral',
  181. 'target_abs_velocity_lateral',
  182. 'courter_angles_speed',
  183. 'target_angles_speed',
  184. 'distance_target',
  185. 'relative_angle_abs_target',
  186. 'relative_angle_abs_courter',
  187. 'relative_orientation_abs_target_wrap',
  188. ]
  189. condition = 'current_song'
  191. start_frame_idx = 2000
  192. end_frame_idx = -2000
  193. # confidence values not calibrated in SLEAP
  194. confidence_threshold = 0.7 #None if sleap
  196. # filtering
  197. required_sampling_rate = 30#Hz
  198. filter_features = False
  199. filter_cutoff = 10#Hz
  200. filter_ord = 4
  201. assert filter_cutoff < required_sampling_rate/2, (
  202. "filter_cutoff should be less than %0.2f Hz"
  203. %(required_sampling_rate/2))
  204. b, a = signal.butter(filter_ord, filter_cutoff/(0.5*required_sampling_rate))
  206. # use normalize_features = True only for wavelet spectrogram approach
  207. normalize_features = True
  208. normalization = "standardize" #"minmax" or None
  209. features_range = None
  211. # time delay works better than wavelet spectrogram for now
  212. time_delay_poses = True
  213. nb_delays = 15
  214. basis_transform = False
  215. basis_function = bases.raised_cosine(0, 12, [0, 12], 10, nb_delays)
  216. wavelet_transform = not time_delay_poses
  217. min_freq = 1
  218. max_freq = 25
  219. nfreqs = 25
  220. freq_spacing = 'dyadic'
  221. freqs, scales_cwt = get_freq_scales(
  222. min_freq, max_freq, nfreqs, required_sampling_rate,
  223. spacing = freq_spacing)
  225. # umap parameters
  226. train_size = 0.1
  227. min_dist = 0.0
  228. n_neighbors = 100
  229. random_state = 42
  230. nan_handling = "remove" # "interpolate" to interpolate nans or "remove" to remove nans
  231. if wavelet_transform:
  232. assert nan_handling == "interpolate", "use nan_handling as interpolate if wavelet_transform is True"
  234. # kde and watershed parameters
  235. bw = 0.5
  236. nb_gridpoints = 256
  237. watershed_threshold = 0.001
  239. results_save_path = f"res/umap_embedding"
  240. if not os.path.exists(results_save_path):
  241. os.makedirs(results_save_path)
  243. # %%
  244. # save umap parameters
  245. umap_params = {}
  246. umap_params = dict(
  247. trial_names=trial_names,
  248. features=features,
  249. condition=condition,
  250. required_sampling_rate=required_sampling_rate,
  251. data_save_path=data_save_path,
  252. start_frame_idx=start_frame_idx,
  253. end_frame_idx=end_frame_idx,
  254. confidence_threshold=confidence_threshold,
  255. filter_features=filter_features,
  256. filter_cutoff=filter_cutoff,
  257. filter_ord=filter_ord,
  258. normalize_features=normalize_features,
  259. normalization=normalization,
  260. features_range=features_range,
  261. #additional_features=additional_features,
  262. time_delay_poses=time_delay_poses,
  263. nb_delays=nb_delays,
  264. basis_transform = basis_transform,
  265. basis_function = basis_function,
  266. wavelet_transform=wavelet_transform,
  267. min_freq=min_freq,
  268. max_freq=max_freq,
  269. nfreqs=nfreqs,
  270. freq_spacing=freq_spacing,
  271. train_size=train_size,
  272. min_dist=min_dist,
  273. n_neighbors=n_neighbors,
  274. random_state=random_state,
  275. nan_handling=nan_handling,
  276. bw=bw,
  277. nb_gridpoints=nb_gridpoints,
  278. watershed_threshold=watershed_threshold,
  279. results_save_path=results_save_path)
  281. # %% [markdown]
  282. # ## Prepare the data
  283. #
  284. # Read the data from `preprocessed_data_path` if available. Otherwise, read the trial-wise CSV files in the `dat` folder. Lowpass-filter the data, remove frames below the confidence threshold, and normalize.
  286. # %%
  287. # get data
  288. all_trial_names = list(np.concatenate(list(trial_names.values())))
  289. X, y = [], []
  290. if os.path.exists(preprocessed_data_path):
  291. data = np.load(preprocessed_data_path, allow_pickle=True)
  292. X = list(data['X'])
  293. y = list(data['y'])
  294. else:
  295. for trial_name in all_trial_names:
  296. print(trial_name)
  297. df_trial_data = pd.read_csv(f'{data_save_path}/{trial_name}.csv')
  298. trial_X = df_trial_data[features].values[start_frame_idx:end_frame_idx]
  299. trial_y = df_trial_data[condition].values[start_frame_idx:end_frame_idx]
  300. # low pass filter data
  301. trial_X = signal.filtfilt(b, a, trial_X, axis=0)
  302. # filter out non confident frames
  303. if confidence_threshold is not None and 'poses_confidence' in df_trial_data:
  304. poses_confidence = df_trial_data['poses_confidence'].values
  305. confident_frames = poses_confidence>confidence_threshold
  306. trial_X[~confident_frames] = np.nan
  307. trial_y[~confident_frames] = np.nan
  308. # normalization
  309. if normalize_features:
  310. if normalization == 'standardize':
  311. trial_X = StandardScaler().fit_transform(trial_X)
  312. elif normalization == 'minmax':
  313. trial_X = MinMaxScaler().fit_transform(trial_X)
  314. X.append(trial_X)
  315. y.append(trial_y)
  316. np.savez(preprocessed_data_path, X=X, y=y)
  318. # %% [markdown]
  319. # ## Data processing
  320. #
  321. # Extract wavelet spectrograms or perform time-delay embedding of input features for each trial. If using wavelet spectrograms, the processed data for each trial is a 2D array of shape (`N_trial`, `len(features)` x `nfreqs`), where `N_trial` is the number of samples for the trial. If using time-delay embedding, the processed data for each trial is a 2D array of shape (`N_trial`, `len(features)` x `nb_delays`). The processed data from all trials are then pooled to form `X_proc`.
  322. #
  324. # %%
  325. # process data
  326. X_proc, y_proc = [], []
  327. list_timestamps = []
  328. for trial_X, trial_y in zip(X, y):
  329. N = trial_X.shape[0]
  330. t0=0
  331. dt=1/required_sampling_rate
  332. timestamps = np.arange(0, N, dtype="float32") * dt + t0
  333. if wavelet_transform:
  334. # Compute wavelet spectrogram
  335. print('Computing wavelet spectrograms ...')
  336. trial_X_cwt = []
  337. for feat_idx in range(trial_X.shape[1]):
  338. [P, t, f] = compute_wavelet(timestamps, trial_X[:, feat_idx],
  339. scales_cwt)
  340. trial_X_cwt.extend((P))
  341. trial_X_cwt = np.array(trial_X_cwt, dtype='float32').T
  342. X_proc.append(trial_X_cwt)
  343. y_proc.append(trial_y)
  344. print('done')
  345. elif time_delay_poses:
  346. # time delay embedding
  347. X_td, y_td = time_delay_embedding(
  348. trial_X, trial_y, nb_delays = nb_delays, multi_features=True,
  349. padding='same', remove_nans=True if nan_handling == "remove" else False)
  350. if basis_transform:
  351. X_td, basis_projection = basis_transformation(
  352. X_td, nb_delays, basis_function, multi_features=True, nb_stim=len(features))
  353. X_proc.append(np.array(X_td, dtype="float32"))
  354. y_proc.append(y_td)
  355. list_timestamps.append(timestamps)
  357. # %% [markdown]
  358. # ## UMAP embedding and spatial clustering using KDE and watershed segmentation
  359. #
  360. # First, the time-delayed inputs `X_proc` are embedded into a two-dimensional manifold `X_embedded`. A kernel density estimation is performed on the two-dimensional space, followed by watershed segmentation, which creates spatial clusters centered at local peaks of the kernel density estimates. From a behavioral perspective, the local peaks in the KDE correspond to stereotyped interactions between the organisms, and regions of low density correspond to transitions between stereotyped interactions. Thus, each spatial cluster is assigned to a particular interaction prototype called a social mode.
  362. # %%
  363. if os.path.exists(f'{results_save_path}/embedding_results.npz'):
  364. # load data
  365. results = np.load(f'{results_save_path}/embedding_results.npz', allow_pickle=True)
  366. X_proc = results['X_proc']
  367. y_proc = results['y_proc']
  368. X_embedded = results['X_embedded']
  369. X_embedded_kde = results['X_embedded_kde']
  370. X_embedded_kde_positions = results['X_embedded_kde_positions']
  371. X_embedded_segments = results['X_embedded_segments']
  372. labels = results['labels']
  373. positions = results['positions']
  374. labels_edge_positions = results['labels_edge_positions']
  375. else:
  376. #%% umap embedding
  377. X_embedded, reducer = umap_embedding(
  378. X_proc, min_dist, n_neighbors,
  379. n_components=2,
  380. train_size=train_size,
  381. random_state=random_state,
  382. )
  383. # kernel density estimation and watershed segmentation
  384. spatial_clustering = SpatialClustering(bw, nb_gridpoints, watershed_threshold)
  385. embedded_features_kde, positions, labels, labels_edge, labels_edge_positions = spatial_clustering.fit(X_embedded)
  386. X_embedded_kde, X_embedded_kde_positions, X_embedded_segments = spatial_clustering.transform(
  387. X_embedded, positions, labels)
  389. # save results
  390. results = {}
  391. results['X_proc'] = X_proc
  392. results['y_proc'] = y_proc
  393. results['timestamps'] = list_timestamps
  394. results['reducer'] = reducer
  395. results['X_embedded'] = X_embedded
  396. results['X_embedded_kde'] = X_embedded_kde
  397. results['X_embedded_kde_positions'] = X_embedded_kde_positions
  398. results['X_embedded_segments'] = X_embedded_segments
  399. results['labels'] = labels
  400. results['positions'] = positions
  401. results['labels_edge_positions'] = labels_edge_positions
  402. results['umap_params'] = umap_params
  403. np.savez(f'{results_save_path}/embedding_results.npz', results)
  405. # %% [markdown]
  406. # ## Meaning of each spatial segment (social modes)
  407. #
  408. # To understand what each spatial segment means, we first plot the time-delay embedded features for each social mode separately for male-female and male-male interactions.
  410. # %%
  411. # time-delay embedded feature values for each segment
  412. X_proc_mean = {}
  413. n_samples_label = {}
  414. for analysis_group in trial_names:
  415. X_proc_mean[analysis_group] = {}
  416. n_samples_label[analysis_group] = {}
  417. for label in np.unique(labels):
  418. if label == 0:
  419. continue
  420. X_proc_mean[analysis_group][label] = []
  421. n_samples_label[analysis_group][label] = 0
  422. for x_proc, x_seg, trial_id in zip(X_proc, X_embedded_segments, all_trial_names):
  423. if trial_id not in trial_names[analysis_group]: continue
  424. X_proc_mean[analysis_group][label].append(
  425. np.sum(x_proc[x_seg==label], axis=0)/\
  426. np.sum(x_seg==label))
  427. n_samples_label[analysis_group][label] += np.sum(x_seg==label)
  428. X_proc_mean[analysis_group][label] = np.array(
  429. X_proc_mean[analysis_group][label])
  430. if time_delay_poses:
  431. X_proc_mean[analysis_group][label] = np.reshape(
  432. X_proc_mean[analysis_group][label],
  433. (X_proc_mean[analysis_group][label].shape[0], -1, nb_delays))
  434. elif wavelet_transform:
  435. X_proc_mean[analysis_group][label] = np.reshape(
  436. X_proc_mean[analysis_group][label],
  437. (X_proc_mean[analysis_group][label].shape[0], -1, nfreqs))
  439. num_clusters = 12
  440. fig, ax = plt.subplots(
  441. len(X_proc_mean),
  442. num_clusters,
  443. sharex=True, sharey=True,
  444. num="cluster_means", figsize=(num_clusters*2, 2*len(X_proc_mean)))
  445. for g, group in enumerate(X_proc_mean):
  446. for i, (label, features_mean) in enumerate(X_proc_mean[group].items()):
  447. ax[g, i].set_title(f"{label} ({n_samples_label[group][label]})")
  448. ax[g, i].matshow(np.nanmean(features_mean, 0), vmin=-1, vmax=1, cmap="bwr")
  449. ax[g, i].set_yticks(np.arange(features_mean.shape[1]))
  450. ax[g, i].set_yticklabels(features)
  451. ax[g, 0].set_ylabel(group)
  452. plt.tight_layout()
  453. plt.show()
  455. # %% [markdown]
  456. # Next, we plot the average values of each input feature within each segment. Based on the speed components of each fly and their relative positioning with respect to each other, we name the modes:
  457. # 1. `Behind and idle`: the courter is behind the target and idle
  458. # 2. `Behind and close`: the courter is behind the target and close
  459. # 3. `Behind and chasing`: the courter is behind the target and chasing
  460. # 4. `Behind and circling`: the courter is behind the target and circling
  461. # 5. `Uninterested`: the courter and target are distant and facing away from each other
  462. # 6. `Front and circling`: the courter is in front of the target and circling
  463. # 7. `Front and close`: the courter is in front of the target and close
  464. # 8. `Front and idle`: the courter is in front of the target and idle
  465. #
  466. # Some segments contain only a few frames or are noisy, which we ignore.
  468. # %%
  469. # feature means for each segment
  470. filter_by_song = False
  471. ignore_states = [0, 1, 5, 9, 12]
  472. state_names = [
  473. 'Behind idle',
  474. 'Behind close',
  475. 'Behind chasing',
  476. 'Behind circling',
  477. 'Uninterested',
  478. 'Front circling',
  479. 'Front close',
  480. 'Front idle'
  481. ]
  482. rows = []
  483. for analysis_group in trial_names:
  484. for trial_name in trial_names[analysis_group]:
  485. trial_idx = trial_names[analysis_group].index(trial_name)
  486. trial_feature_values = StandardScaler().fit_transform(X[trial_idx])
  487. trial_segments = X_embedded_segments[trial_idx]
  488. for state in np.unique(trial_segments):
  489. if state in ignore_states: continue
  490. for f, feat in enumerate(features):
  491. trial_feat_segment_mean = np.nanmean(
  492. trial_feature_values[trial_segments==state, f])
  493. rows.append(
  494. {"trial_name": trial_name,
  495. "analysis_group": analysis_group,
  496. "state": state,
  497. "feature": feat,
  498. "zscore": trial_feat_segment_mean}
  499. )
  500. df_mean_feature_values_segments = pd.DataFrame(rows)
  501. # plots
  502. figname = 'cluster_feature_means'
  503. fig, ax = plt.subplots(
  504. 1, len(np.unique(df_mean_feature_values_segments.state)),
  505. figsize=(3*(len(np.unique(df_mean_feature_values_segments.state))),
  506. len(np.unique(df_mean_feature_values_segments.feature))),
  507. sharey=True,
  508. num=figname)
  509. for s, state in enumerate(np.unique(df_mean_feature_values_segments.state)):
  510. if state in ignore_states: continue
  511. df_mean_feature_values = df_mean_feature_values_segments[
  512. df_mean_feature_values_segments.state==state]
  513. sns.barplot(data=df_mean_feature_values,
  514. y="feature",
  515. x="zscore",
  516. ax=ax[s],
  517. dodge=True,
  518. #palette=["k"],
  519. orient='h')
  520. #sns.pointplot(data=df_mean_feature_values,
  521. # y="feature",
  522. # x="zscore",
  523. # dodge=0.5,
  524. # ax=ax[s],
  525. # orient='h',
  526. # join=False, palette=['k'])
  527. ax[s].set_title(f"{state_names[s]}")
  528. ax[s].set_yticklabels(features)
  529. ax[s].legend([])
  530. plt.tight_layout()
  531. plt.show()
  534. # %% [markdown]
  535. # Next, we visualize the feature values directly on the social UMAP embedding.
  537. # %%
  538. # color code embedding by feature values
  539. binned = True
  540. filter_song = False
  541. groupwise=False
  542. figname = "embedding_feat_values"
  543. if binned: figname = figname + "_binned"
  544. if groupwise: figname = figname + "_groupwise"
  545. if groupwise:
  546. feat_values_binned = {}
  547. fig, ax = plt.subplots(
  548. len(trial_names), len(features), sharex=True, sharey=True,
  549. num='embedding_feature_values'+("_binned" if binned else "_"),
  550. figsize=(len(features)*2.5, 6))
  551. for g, analysis_group in enumerate(trial_names):
  552. feat_values_binned[analysis_group] = {}
  553. X_embedded_concat_group = np.concatenate(
  554. [X_embedded[i] for i in range(len(X_embedded))
  555. if all_trial_names[i] in trial_names[analysis_group]], axis=0)
  556. X_raw_concat_group = np.concatenate(
  557. [X[i] for i in range(len(X))
  558. if all_trial_names[i] in trial_names[analysis_group]], axis=0)
  559. X_raw_concat_group = StandardScaler().fit_transform(X_raw_concat_group)
  560. for i, feature in enumerate(features):
  561. if binned:
  562. ret = binned_statistic_2d(
  563. X_embedded_concat_group[:, 1],
  564. X_embedded_concat_group[:, 0],
  565. X_raw_concat_group[:, i],
  566. np.nanmean, bins=128,
  567. range=[[0, 18], [0, 18]])
  568. feat_values_binned[analysis_group][feature] = ret
  569. im = ax[g, i].imshow(
  570. ret.statistic,
  571. extent=[ret.x_edge[0], ret.x_edge[-1],
  572. ret.y_edge[0], ret.y_edge[-1]],
  573. cmap='bwr', zorder=1, vmin=-1, vmax=1,)
  574. else:
  575. ax[g, i].scatter(
  576. *X_embedded_concat_group[::10].T, s=0.1, vmin=-1, vmax=1,
  577. alpha=0.25, c=X_raw_concat_group[::10, i], cmap='bwr', zorder=1)
  578. ax[g, i].scatter(*labels_edge_positions.T, s=0.2, c='k', zorder=2)
  579. ax[g, i].set_title(feature)
  580. ax[g, i].set_xlim(0, 18)
  581. ax[g, i].set_ylim(0, 18)
  582. ax[g, 0].set_ylabel(analysis_group)
  583. ax[g, i].axis("off")
  584. fig.colorbar(
  585. im, ax=ax[-1, -1], shrink=0.25, aspect=10, ticks=[-1, 0, 1])
  586. # both groups together
  587. else:
  588. fig, ax = plt.subplots(
  589. 1, len(features), sharex=True, sharey=True,
  590. num='embedding_feature_values'+("_binned" if binned else "_"),
  591. figsize=(len(features)*2.5, 3))
  592. feat_values_binned = {}
  593. if not filter_song:
  594. X_embedded_concat_group = np.concatenate(
  595. [X_embedded[i] for i in range(len(X_embedded))], axis=0)
  596. X_raw_concat_group = np.concatenate(
  597. [X[i] for i in range(len(X))], axis=0)
  598. else:
  599. X_embedded_concat_group = np.concatenate(
  600. [X_embedded[i][y[i].ravel()!=0]
  601. for i in range(len(X_embedded))], axis=0)
  602. X_raw_concat_group = np.concatenate(
  603. [X[i][y[i].ravel()!=0] for i in range(len(X))], axis=0)
  604. X_raw_concat_group = StandardScaler().fit_transform(X_raw_concat_group)
  605. for i, feature in enumerate(features):
  606. if binned:
  607. ret = binned_statistic_2d(
  608. X_embedded_concat_group[:, 1],
  609. X_embedded_concat_group[:, 0],
  610. X_raw_concat_group[:, i],
  611. np.nanmean, bins=128,
  612. range=[[-18, 18], [-18, 18]])
  613. feat_values_binned[feature] = ret
  614. im = ax[i].imshow(
  615. ret.statistic,
  616. extent=[ret.x_edge[0], ret.x_edge[-1],
  617. ret.y_edge[0], ret.y_edge[-1]],
  618. cmap='bwr', zorder=1, vmin=-1, vmax=1,)
  619. else:
  620. im = ax[i].scatter(
  621. *X_embedded_concat_group[::10].T, s=0.1, vmin=-1, vmax=1,
  622. alpha=0.25, c=X_raw_concat_group[::10, i], cmap='bwr', zorder=1)
  623. ax[i].scatter(*labels_edge_positions.T, s=0.2, c='k', zorder=2)
  624. ax[i].set_title(feature, fontsize=10)
  625. ax[i].set_xlim(0, 20)
  626. ax[i].set_ylim(0, 20)
  627. ax[i].axis("off")
  628. fig.colorbar(im, ax=ax[-1], shrink=0.25, aspect=10, ticks=[-1, 0, 1])
  629. plt.tight_layout()
  630. plt.show()
  632. # %% [markdown]
  633. # Finally, we plot the courting male's position around the partner during each social mode. We normalize the distance between flies by the length of the target and limit the visualization to three fly lengths.
  635. # %%
  636. if os.path.exists(f'{results_save_path}/courter_positions.npz'):
  637. with np.load(f'{results_save_path}/courter_positions.npz', allow_pickle=True) as data:
  638. dict_courter_positions = data['arr_0'].item()
  639. else:
  640. dict_courter_positions = {}
  641. for trial_idx, trial_name in enumerate(all_trial_names):
  642. df_trial_data = pd.read_csv(f'{data_save_path}/{trial_name}.csv')
  643. courter_position = df_trial_data[['relative_angle_courter', 'distance_target']].values
  644. # normalize distance by target length
  645. courter_position[:, 1] = courter_position[:, 1]/df_trial_data["target_length"].values
  646. trial_segments = X_embedded_segments[trial_idx]
  647. for state in np.unique(trial_segments):
  648. if state in ignore_states: continue
  649. else:
  650. if state not in dict_courter_positions:
  651. dict_courter_positions[state] = []
  652. dict_courter_positions[state].extend(courter_position[:len(trial_segments)][trial_segments==state])
  653. np.savez(f'{results_save_path}/courter_positions.npz', dict_courter_positions)
  655. rbins = np.linspace(0, 3, 31)
  656. abins = np.linspace(-np.pi, np.pi, 73)
  657. courter_position_hist = {}
  658. for state in dict_courter_positions:
  659. dict_courter_positions[state]=np.array(dict_courter_positions[state])
  660. azimut = np.array(dict_courter_positions[state][:, 0])
  661. azimut_rad = np.deg2rad(azimut)
  662. radius = np.array(dict_courter_positions[state][:, 1])
  663. #calculate histogram
  664. hist, _, _ = np.histogram2d(azimut_rad, radius, bins=(abins, rbins),
  665. density=True)
  666. courter_position_hist[state] = hist
  668. # plot
  669. A, R = np.meshgrid(abins, rbins)
  670. figname = 'courter_position_states'
  671. fig, ax = plt.subplots(
  672. 1, len(courter_position_hist), subplot_kw=dict(projection="polar"),
  673. figsize=(3*len(courter_position_hist), 3), num=figname)
  674. for s, (segment, hist) in enumerate(courter_position_hist.items()):
  675. pc = ax[s].pcolormesh(
  676. A, R, (hist).T,
  677. cmap="Reds",
  678. vmin=0, vmax=1)
  679. ax[s].set_title(state_names[s], fontsize=15)
  680. ax[s].set_xticks([-np.pi*(3/4), -np.pi/2, -np.pi/4,
  681. 0, np.pi/4, np.pi/2, np.pi*(3/4), np.pi])
  682. ax[s].set_xlim(-np.pi, np.pi)
  683. ax[s].set_theta_zero_location("N")
  684. plt.subplots_adjust(bottom=0.4, right=0.9, top=0.6)
  685. cax = plt.axes([0.85, 0.1, 0.01, 0.2])
  686. plt.colorbar(pc, cax=cax)
  687. plt.tight_layout()
  688. plt.show()
  691. # %% [markdown]
  692. # ## Comparison between male-female and male-male interactions
  693. #
  694. # Once we have the social UMAP embedding, we can compare them across different experimental conditions to understand differences in their interactions. Here we plot the UMAPs for male-female and male-male interactions and visualize their differences.
  696. # %%
  697. # comparison between groups
  698. kde_mean_group = {}
  699. X_embedded_kde_group = {}
  700. groups = list(trial_names.keys())
  701. for group in groups:
  702. X_embedded_kde_group[group] = []
  703. for trial_name in trial_names[group]:
  704. trial_idx = all_trial_names.index(trial_name)
  705. X_embedded_kde_group[group].append(X_embedded_kde[trial_idx])
  706. group_kde_mean = np.mean(X_embedded_kde_group[group], 0)
  707. group_kde_mean[group_kde_mean<1e-3]=0
  708. kde_mean_group[group] = group_kde_mean
  710. fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(15, 5))
  711. xmin = positions[:,0].min()
  712. xmax = positions[:,0].max()
  713. ymin = positions[:,1].min()
  714. ymax = positions[:,1].max()
  715. for g, group in enumerate(kde_mean_group):
  716. vmin = 0 #np.round(np.nanpercentile(kde_mean_group[group], 2.5), 2)
  717. vmax = 0.03 #np.round(np.nanpercentile(kde_mean_group[group], 97.5), 2)
  718. kde_plot = ax[g].imshow(
  719. kde_mean_group[group],
  720. origin='lower',
  721. extent=[xmin, xmax, ymin, ymax],
  722. vmin=vmin, vmax=vmax,
  723. cmap="Reds")
  724. ax[g].scatter(*labels_edge_positions.T, s=0.2, c='k')
  725. ax[g].set_title(group)
  726. ax[g].axis('off')
  727. ax[g].set_xlim(0, 20)
  728. ax[g].set_ylim(0, 20)
  729. fig.colorbar(kde_plot, ax=ax[g], shrink=0.25, aspect=10, ticks=[vmin, vmax])
  730. # difference
  731. statistical_test = True
  732. X_embedded_kde_diff = kde_mean_group[groups[0]] - kde_mean_group[groups[1]]
  733. vmin=-0.03
  734. vmax=0.03
  735. if statistical_test:
  736. statistic, p_value = mannwhitneyu(X_embedded_kde_group[groups[0]], X_embedded_kde_group[groups[1]])
  737. X_embedded_kde_diff = X_embedded_kde_diff * (p_value<0.05)
  738. diff_lim = np.max(np.abs(X_embedded_kde_diff))
  739. kde_diff = ax[2].imshow(
  740. X_embedded_kde_diff,
  741. cmap='bwr', origin='lower',
  742. extent=[xmin, xmax, ymin, ymax],
  743. vmin=vmin, vmax=vmax)
  744. ax[2].scatter(*labels_edge_positions.T, s=0.2, c='k')
  745. ax[2].set_title("difference")
  746. ax[2].set_xlim(0, 20)
  747. ax[2].set_ylim(0, 20)
  748. ax[2].axis('off')
  749. fig.colorbar(kde_diff, ax=ax[2], shrink=0.25, aspect=10, ticks=[vmin, 0, vmax])
  750. plt.tight_layout()
  751. plt.show()
  753. plt.tight_layout()
  754. plt.show()
  756. # %% [markdown]
  757. # Next, we quantify the differences in social modes during male-female and male-male interactions.
  759. # %%
  760. rows = []
  761. for analysis_group in trial_names:
  762. for trial_name in trial_names[analysis_group]:
  763. trial_idx = all_trial_names.index(trial_name)
  764. trial_segments = X_embedded_segments[trial_idx]
  765. for s, state in enumerate(np.unique(trial_segments)):
  766. if s in ignore_states: continue
  767. state_ratio = np.sum(trial_segments==state)/len(trial_segments)
  768. rows.append(
  769. {"trial_name": trial_name,
  770. "analysis_group": analysis_group,
  771. "state": s,
  772. "ratio": state_ratio}
  773. )
  775. df_states_time_spent = pd.DataFrame(rows)
  776. fig, ax = plt.subplots(1, 1, figsize=(2*len(np.unique(df_states_time_spent.state)), 4))
  777. sns.barplot(data=df_states_time_spent, x="state", y="ratio", hue="analysis_group")
  778. sns.stripplot(data=df_states_time_spent, x="state", y="ratio", hue="analysis_group", dodge=True, alpha=0.25, palette=["k"])
  779. sns.despine()
  780. ax.set_ylim(0, 1)
  781. ax.set_xticklabels(state_names)
  782. plt.tight_layout()
  783. plt.show()
  786. # %% [markdown]
  787. # ## Comparisons between different song contexts
  788. #
  789. # To compare interactions in different singing conditions, we condition our UMAPs on different song contexts and compare song versus silence and pulse versus sine.
  791. # %%
  792. unique_segments = np.unique(np.concatenate(X_embedded_segments))
  793. song_conditioned_kde = {}
  794. rows = []
  795. for analysis_group in trial_names:
  796. song_conditioned_kde[analysis_group] = {}
  797. for song_type in ['song', 'silence', 'pulse', 'sine']:
  798. song_conditioned_kde[analysis_group][song_type] = []
  799. for datename in trial_names[analysis_group]:
  800. trial_idx = all_trial_names.index(datename)
  801. X_embedded_trial = X_embedded[trial_idx][nb_delays:]
  802. X_embedded_segments_trail = X_embedded_segments[trial_idx][nb_delays:]
  803. y_trial = y_proc[trial_idx].ravel()
  804. # song
  805. if song_type == 'song':
  806. trial_song_embedding = X_embedded_trial[(y_trial!=0)]
  807. trial_song_segments = X_embedded_segments_trail[y_trial!=0]
  808. elif song_type == 'silence':
  809. trial_song_embedding = X_embedded_trial[(y_trial==0)]
  810. trial_song_segments = X_embedded_segments_trail[y_trial==0]
  811. elif song_type == 'pulse':
  812. trial_song_embedding = X_embedded_trial[(y_trial==1)]
  813. trial_song_segments = X_embedded_segments_trail[y_trial==1]
  814. elif song_type == 'sine':
  815. trial_song_embedding = X_embedded_trial[(y_trial==2)]
  816. trial_song_segments = X_embedded_segments_trail[y_trial==2]
  817. _, trial_song_kde = kernel_density_estimation(
  818. trial_song_embedding, bw, positions)
  819. song_conditioned_kde[analysis_group][song_type].append(trial_song_kde)
  820. for state in unique_segments:
  821. if state not in ignore_states:
  822. song_state_time_trial = (
  823. np.sum(trial_song_segments==state) / len(trial_song_segments)
  824. )
  825. rows.append(
  826. {"analysis_group": analysis_group,
  827. "song": song_type,
  828. "trial": datename,
  829. "state": state,
  830. "time_spent": np.sum(
  831. trial_song_segments==state),
  832. "time_spent_ratio": song_state_time_trial}
  833. )
  835. song_conditioned_kde[analysis_group][song_type] = np.array(
  836. song_conditioned_kde[analysis_group][song_type])
  838. df_song_states_time = pd.DataFrame(rows)
  841. # %%
  842. # song vs silence
  843. rows = []
  844. for analysis_group in np.unique(df_song_states_time.analysis_group):
  845. for trial in np.unique(df_song_states_time.trial):
  846. for state in np.unique(df_song_states_time.state):
  847. trial_song_states_time = df_song_states_time[
  848. (df_song_states_time.analysis_group==analysis_group)&
  849. (df_song_states_time.trial==trial)&
  850. (df_song_states_time.state==state)
  851. ]
  852. if not len(trial_song_states_time): continue
  853. song_silence_diff = trial_song_states_time[
  854. (trial_song_states_time.song=="song")].time_spent_ratio.values[0] - trial_song_states_time[
  855. (trial_song_states_time.song=="silence")].time_spent_ratio.values[0]
  856. pulse_sine_diff = trial_song_states_time[
  857. (trial_song_states_time.song=="pulse")].time_spent_ratio.values[0] - trial_song_states_time[
  858. (trial_song_states_time.song=="sine")].time_spent_ratio.values[0]
  859. rows.append(
  860. {"analysis_group": analysis_group,
  861. "trial": trial,
  862. "state": state,
  863. "p(song)-p(silence)": song_silence_diff,
  864. "p(pulse)-p(sine)": pulse_sine_diff}
  865. )
  866. df_song_diff = pd.DataFrame(rows)
  867. fig, ax = plt.subplots(1, 3, figsize=(15, 3), gridspec_kw={'width_ratios': [1, 1, 3]})
  868. mean_kde_mf_diff = np.mean(song_conditioned_kde['male-female']['song'], 0) - \
  869. np.mean(song_conditioned_kde['male-female']['silence'], 0)
  870. mean_kde_mm_diff = np.mean(song_conditioned_kde['male-male']['song'], 0) - \
  871. np.mean(song_conditioned_kde['male-male']['silence'], 0)
  872. mf_kde = ax[0].matshow(mean_kde_mf_diff,
  873. cmap='bwr', origin='lower',
  874. extent=[xmin, xmax, ymin, ymax],
  875. vmin=-0.03, vmax=0.03)
  876. ax[0].set_xlim(0, 20)
  877. ax[0].set_ylim(0, 20)
  878. ax[0].set_ylabel("song - silence")
  879. ax[0].set_title("female-directed")
  880. ax[0].axis("off")
  881. ax[0].scatter(*labels_edge_positions.T, s=0.2, c='k')
  882. fig.colorbar(mf_kde, ax=ax[0], shrink=0.25, aspect=10, ticks=[-0.03, 0, 0.03])
  883. mm_kde = ax[1].matshow(mean_kde_mm_diff,
  884. cmap='bwr', origin='lower',
  885. extent=[xmin, xmax, ymin, ymax],
  886. vmin=-0.03, vmax=0.03)
  887. ax[1].set_xlim(0, 20)
  888. ax[1].set_ylim(0, 20)
  889. ax[1].axis("off")
  890. ax[1].set_title("male-directed")
  891. ax[1].scatter(*labels_edge_positions.T, s=0.2, c='k')
  892. #fig.colorbar(mm_kde, ax=ax[1], shrink=0.25, aspect=10, ticks=[-0.03, 0, 0.03])
  894. sns.barplot(data=df_song_diff, x="state", y="p(song)-p(silence)", hue="analysis_group", ax=ax[2])
  895. sns.stripplot(data=df_song_diff, x="state", y="p(song)-p(silence)", hue="analysis_group", ax=ax[2], alpha=0.25, dodge=True, palette=["k"])
  896. ax[2].set_xticklabels(state_names, rotation=45)
  897. plt.suptitle("song - silence")
  898. #plt.tight_layout()
  901. # %%
  902. # pulse vs sine
  903. fig, ax = plt.subplots(1, 3, figsize=(15, 3), gridspec_kw={'width_ratios': [1, 1, 3]})
  904. mean_kde_mf_diff = np.mean(song_conditioned_kde['male-female']['pulse'], 0) - \
  905. np.mean(song_conditioned_kde['male-female']['sine'], 0)
  906. mean_kde_mm_diff = np.mean(song_conditioned_kde['male-male']['pulse'], 0) - \
  907. np.mean(song_conditioned_kde['male-male']['sine'], 0)
  908. mf_kde = ax[0].matshow(mean_kde_mf_diff,
  909. cmap='bwr', origin='lower',
  910. extent=[xmin, xmax, ymin, ymax],
  911. vmin=-0.03, vmax=0.03)
  912. ax[0].set_xlim(0, 20)
  913. ax[0].set_ylim(0, 20)
  914. ax[0].set_title("female-directed")
  915. ax[0].axis("off")
  916. ax[0].scatter(*labels_edge_positions.T, s=0.2, c='k')
  917. fig.colorbar(mf_kde, ax=ax[0], shrink=0.25, aspect=10, ticks=[-0.03, 0, 0.03])
  918. mm_kde = ax[1].matshow(mean_kde_mm_diff,
  919. cmap='bwr', origin='lower',
  920. extent=[xmin, xmax, ymin, ymax],
  921. vmin=-0.03, vmax=0.03)
  922. ax[1].set_xlim(0, 20)
  923. ax[1].set_ylim(0, 20)
  924. ax[1].axis("off")
  925. ax[1].set_title("male-directed")
  926. ax[1].scatter(*labels_edge_positions.T, s=0.2, c='k')
  927. #fig.colorbar(mm_kde, ax=ax[1], shrink=0.25, aspect=10, ticks=[-0.03, 0, 0.03])
  929. sns.barplot(data=df_song_diff, x="state", y="p(pulse)-p(sine)", hue="analysis_group", ax=ax[2])
  930. sns.stripplot(data=df_song_diff, x="state", y="p(pulse)-p(sine)", hue="analysis_group", ax=ax[2], alpha=0.25, dodge=True, palette=["k"])
  931. ax[2].set_xticklabels(state_names, rotation=45)
  932. plt.suptitle("pulse - sine")
  933. #plt.tight_layout()
  935. # %% [markdown]
  936. # ## Hyperparameter tuning
  937. #
  938. # We used three important hyperparameters: two for fitting the UMAP, `n_neighbors` and `min_dist`, and one for time-delay embedding of features, `nb_delays`.
  939. #
  940. # The `n_neighbors` parameter specifies the size of the local neighborhood UMAP looks at when fitting the manifold. This parameter controls how UMAP balances local versus global structure in the data. Low values make the UMAP fit the local structure well, whereas large values focus on global structure but can lose fine details.
  941. #
  942. # `min_dist` specifies the minimum distance between the embedded points in the low-dimensional manifold. Thus, low values result in points being embedded densely together.
  943. #
  944. # `nb_delays` specifies the history information to be included when embedding the features into low dimensions. Small values focus on immediate history.
  945. #
  946. # To find out the optimal values of these hyperparameters for our data, we performed a hyperparameter optimization. This was done by fitting the UMAP with different hyperparameter combinations and validating performance by the reconstruction error. From the fitted UMAP, an `inverse_transform` is applied on an embedded validation set (not used for fitting), and the mean squared error between the original validation data and reconstructed validation data is quantified.
  948. # %%
  949. #UMAP hyperparameter tuning
  950. n_neighbors_grid=[10, 50, 100, 200]
  951. min_dist_grid=[0.0, 0.1, 0.2, 0.5]
  952. nb_delays_grid=[15, 30, 60, 120]#s
  953. fit_results = umap_multiple_fits(
  954. X,
  955. n_neighbors_grid=n_neighbors_grid,
  956. min_dist_grid=min_dist_grid,
  957. nb_delays_grid=nb_delays_grid)
  959. # plots
  960. val_score_matrices = {}
  961. for nb_delays in nb_delays_grid:
  962. val_score_matrices[nb_delays] = np.zeros((len(min_dist_grid), len(n_neighbors_grid)))
  963. for m, min_dist in enumerate(min_dist_grid):
  964. for n, n_neighbors in enumerate(n_neighbors_grid):
  965. val_score_matrices[nb_delays][m, n] = fit_results[
  966. f'delay: {nb_delays}, '+\
  967. f'n_components: 2, ' +\
  968. f'n_neighbors: {n_neighbors}, ' +\
  969. f'min_dist: {min_dist}, '+\
  970. f'metric: euclidean'
  971. ]['val_score']
  973. fig, ax = plt.subplots(1, len(nb_delays_grid), figsize=(len(nb_delays_grid)*3, 3))
  974. for i, nb_delays in enumerate(nb_delays_grid):
  975. im = ax[i].matshow(val_score_matrices[nb_delays])
  976. ax[i].set_xlabel("n_neighbors")
  977. ax[i].set_xticklabels(n_neighbors_grid)
  978. ax[i].set_ylabel("min_dist")
  979. ax[i].set_yticklabels(min_dist_grid)
  980. fig.colorbar(im, ax=ax[i], shrink=0.25, aspect=10)
  981. plt.tight_layout()
  982. plt.show()
  984. # %%

social_behavior_umap.ipynb at commit db64ff6, no license · at the source

Overview

Authors: Sarath Ravindran Nair1, Adrián Palacios-Muñoz1,2,3, Sage Martineau2,4, Malak Nasr2,5,6, Jan Clemens1,2,7
ORCID iDs: Adrián Palacios-Muñoz, Jan Clemens
  1. ENI-G, a Joint Initiative of the University Medical Center Göttingen and the Max Planck Institute for Multidisciplinary Sciences, Göttingen, Germany
  2. IMPRS Neuroscience, Göttingen, Germany
  3. Present Address: Charité Universitätsmedizin Berlin, Berlin, Germany
  4. Institute of Computer Science, University of Göttingen, Göttingen, Germany
  5. CERVO Brain Research Centre, Québec City, QC Canada
  6. Faculty of Medicine, Université Laval, Québec City, QC Canada
  7. Department of Neuroscience, Faculty VI, University of Oldenburg, Oldenburg, Germany
Institutions: Universitätsmedizin Göttingen (Germany); European Neuroscience Institute Göttingen (Germany); Max Planck Institute for Multidisciplinary Sciences (Germany); University of Göttingen (Germany); Charité - Universitätsmedizin Berlin (Germany); Université Laval (Canada); Carl von Ossietzky Universität Oldenburg (Germany)
Journal: Nature communications, volume 17, issue 1, article 4026
Dates: received 21 July 2025; accepted 2 April 2026; published online 4 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-72057-9 · PMID 42082469 · PMCID PMC13139495 · OpenAlex W4410238817
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.05.08.652884
Status: code verified
Categories: drosophila (organism), systems (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Spectral & time-frequency, Preprocessing, Evoked potentials, fMRI & imaging
Keywords: Social behaviour, Decision, Sexual behaviour
MeSH: Drosophila melanogaster*, Sexual Behavior, Animal*, Social Behavior*, Animals, Brain, Courtship, Female, Male, Optogenetics, Vocalization, Animal (* major topic)
Topic: Neuroendocrine regulation and behavior (Social Psychology, Psychology), according to OpenAlex
Funding: European Research Council (851210); Deutsche Forschungsgemeinschaft (German Research Foundation) (329518246)
Citations: not cited yet (Europe PMC); 138 references in the paper

Abstract

How the brain enables individuals to adapt behavior to their partner is key to understanding social exchange. For example, courtship behavior involves sensorimotor processing of signals that can result in behavioral dialog between partners, such as stereotyped movements and singing. The courtship behavior of Drosophila melanogaster males with their partners, which are usually female but can also be male, involves singing. To investigate how behavioral feedback and sensorimotor processing contribute to flexible social interactions, we compared the courtship behavior and singing of male D. melanogaster towards males and females. Quantitative analysis of their interactions revealed that while underlying courtship and song rules are unaffected by the sex of the partner, the behavioral dynamics and song sequences differ by partner sex. This divergence stems from sex-specific behavioral feedback: females decelerate to song, while males orient towards the singer. Moreover, optogenetic manipulations reveal that the partners’ responses are driven by sex-specific neural circuits that link song detection with arousal and social decisions. Our findings demonstrate that flexible social behaviors can arise from fixed sensorimotor rules through a context-dependent selection facilitated by the partner’s behavioral feedback. More broadly, our results reveal compositionality as a key mechanism for achieving behavioral flexibility during complex social interactions such as courtship.

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

Repositories

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

janclemenslab/socialUMAP

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: db64ff6a2513c17db493d9c7d931daea4bce087b, 26 January 2026
Languages: Python (1), Jupyter (1)
Size: 10 files, 2 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (requirements.txt), 1 notebook
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (2 files), NumPy (2 files), pandas (2 files), scikit-learn (2 files), SciPy (2 files), PyWavelets (1 file), scikit-image (1 file), seaborn (1 file), UMAP (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
3 files
At the source: github.com/janclemenslab/socialUMAP

jgraving/DeepPoseKit

License: Apache-2.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: cecdb0c8c364ea049a3b705275ae71a2f366d4da, 9 June 2020
Languages: Python (52), Jupyter (7)
Size: 66 files, 59 scripts
Software Heritage: not archived
Found in: “Software and algorithms used”
Holds: README, license file, environment (setup.py), 7 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (35 files), TensorFlow (30 files), Keras (26 files), OpenCV (11 files), Matplotlib (9 files), h5py (8 files), pandas (5 files), scikit-learn (1 file), SciPy (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
61 files
At the source: github.com/jgraving/DeepPoseKit

janclemenslab/das

License: Apache-2.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 9a7b5b5ac850a7f85d05bd0f88477ea2991f46d2, 15 September 2026
Languages: Python (48), Jupyter (11)
Size: 125 files, 59 scripts
Software Heritage: archived
Found in: “Software and algorithms used”
Holds: README, license file, environment (pyproject.toml), tests, continuous integration, documentation, 11 notebooks
Not found: CITATION.cff
Tools: NumPy (38 files), Keras (20 files), Matplotlib (10 files), SciPy (9 files), h5py (7 files), pandas (6 files), scikit-learn (6 files), UMAP (3 files), PyTorch (2 files), xarray (2 files), seaborn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
61 files
At the source: github.com/janclemenslab/das

janclemenslab/glm_utils

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 65b048a1d9f2aa0a8f21be2d78ab154ecd05f241, 8 June 2021
Languages: Jupyter (7), Python (6)
Size: 19 files, 13 scripts
Software Heritage: archived
Found in: “Software and algorithms used”
Holds: README, license file, environment (pyproject.toml, setup.py), 7 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (12 files), scikit-learn (8 files), Matplotlib (6 files), SciPy (6 files), h5py (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
15 files
At the source: github.com/janclemenslab/glm_utils

lindermanlab/ssm

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: eb6c8aa33e5311d3564075807dec340759dd8081, 9 May 2025
Languages: Python (51), Jupyter (14)
Size: 93 files, 65 scripts
Software Heritage: archived
Found in: “Software and algorithms used”
Holds: README, license file, CITATION.cff, environment (setup.cfg, setup.py), tests, continuous integration, documentation, 14 notebooks
Tools: Matplotlib (36 files), seaborn (25 files), SciPy (17 files), NumPy (11 files), scikit-learn (11 files), Numba (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
67 files
At the source: github.com/lindermanlab/ssm/

lmcinnes/umap

License: BSD-3-Clause
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 1180b785023ff8e8eb071de56c581d1dc8bec04e, 24 September 2026
Languages: Python (48), Jupyter (11), Shell (3)
Size: 319 files, 62 scripts
Software Heritage: archived
Found in: “Software and algorithms used”
Holds: README, license file, environment (pyproject.toml, setup.py), tests, continuous integration, documentation, 11 notebooks
Not found: CITATION.cff
Tools: UMAP (50 files), NumPy (44 files), scikit-learn (38 files), Matplotlib (20 files), Numba (20 files), SciPy (15 files), TensorFlow (10 files), Keras (9 files), pandas (6 files), seaborn (6 files), h5py (1 file), Pillow (1 file), PyTorch (1 file), scikit-image (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
64 files
At the source: github.com/lmcinnes/umap

tommyod/KDEpy

License: BSD-3-Clause
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: cdeca713299f17ab0a55365c761c1999dc0421a0, 6 May 2026
Languages: Python (26), MATLAB (9), Jupyter (6)
Size: 76 files, 41 scripts
Software Heritage: not archived
Found in: “Software and algorithms used”
Holds: README, license file, environment (pyproject.toml, setup.py), tests, continuous integration, documentation, 6 notebooks
Not found: CITATION.cff
Tools: NumPy (30 files), SciPy (13 files), Matplotlib (12 files), scikit-learn (3 files), statsmodels (3 files), Numba (1 file), pandas (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
43 files
At the source: github.com/tommyod/KDEpy

trevismd/statannotations

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 3f020ae631ca88a091b6ee3e9a9fd32158920879, 22 June 2025
Languages: Python (37), JavaScript (9), Shell (2), Jupyter (1)
Size: 119 files, 49 scripts
Software Heritage: archived
Found in: “Software and algorithms used”
Holds: README, license file, environment (requirements.txt, setup.py, docs/requirements.in, docs/requirements.txt), tests, continuous integration, documentation, 1 notebook
Not found: CITATION.cff
Tools: statannotations (30 files), pandas (14 files), seaborn (14 files), NumPy (13 files), Matplotlib (7 files), SciPy (5 files), statsmodels (2 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
51 files
At the source: github.com/trevismd/statannotations

pywavelets.readthedocs.io

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: the link answers
Software Heritage: not checked
Found in: “Software and algorithms used”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
At the source: pywavelets.readthedocs.io

Code availability

All data analyses were performed using the software listed in Table 2. Code for generating the social maps is deposited at https://github.com/janclemenslab/socialUMAP.

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:

  • 9 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 350 scripts, each with its path and the digest of its content;
  • 11 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

All data supporting the findings of this study are available within the paper, its Supplementary Information or a public repository. Source data are provided with this paper as a Source Data file. Raw experimental data generated in this study have been deposited in the Göttingen Research Online database (10.25625/PTKYTD). 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, 28 September 2026: the first record

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

Cite

This paper

Ravindran Nair, S., Palacios-Muñoz, A., Martineau, S., Nasr, M., & Clemens, J. (2026). Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior. Nature communications, 17(1), 4026. https://doi.org/10.1038/s41467-026-72057-9

BibTeX

@article{ravindrannair2026sex,
author = {Ravindran Nair, Sarath and Palacios-Muñoz, Adrián and Martineau, Sage and Nasr, Malak and Clemens, Jan},
title = {{Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior}},
journal = {Nature communications},
year = {2026},
month = may,
volume = {17},
number = {1},
pages = {4026},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-72057-9},
url = {https://doi.org/10.1038/s41467-026-72057-9},
pmid = {42082469},
pmcid = {PMC13139495}
}

RIS

TY - JOUR
AU - Ravindran Nair, Sarath
AU - Palacios-Muñoz, Adrián
AU - Martineau, Sage
AU - Nasr, Malak
AU - Clemens, Jan
TI - Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/05/04
VL - 17
IS - 1
SP - 4026
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-72057-9
UR - https://doi.org/10.1038/s41467-026-72057-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-72057-9",
"type": "article-journal",
"title": "Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior",
"container-title": "Nature communications",
"author": [
{
"family": "Ravindran Nair",
"given": "Sarath"
},
{
"family": "Palacios-Muñoz",
"given": "Adrián"
},
{
"family": "Martineau",
"given": "Sage"
},
{
"family": "Nasr",
"given": "Malak"
},
{
"family": "Clemens",
"given": "Jan"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "4026",
"DOI": "10.1038/s41467-026-72057-9",
"PMID": "42082469",
"PMCID": "PMC13139495",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-72057-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
4
]
]
}
}

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.3389/fnsys.2026.1822122 [code]
Convergence-divergence circuits for multimodal integration of innate and learned opponent valences.
Journal: Frontiers in systems neuroscience
In common: Keras, UMAP, Numba, 11 other tools, systems, 13 references
[2] doi:10.1038/s41586-026-10735-w [code]
Distributed control circuits across a brain-and-cord connectome.
Journal: Nature
In common: UMAP, scikit-image, PyTorch, 6 other tools, drosophila, 9 references
[3] doi:10.1093/g3journal/jkag037 [code]
Visual recognition of the anteroposterior female body axis drives spatial elements of male courtship in Drosophila.
Journal: G3 (Bethesda, Md.)
In common: scikit-image, h5py, scikit-learn, 4 other tools, drosophila, 8 references
[4] doi:10.1038/s41592-026-03057-2 [code]
CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: PyWavelets, Keras, UMAP, 12 other tools
[5] doi: [code]
Real-time closed-loop feedback system for mouse mesoscale cortical signal and movement control
Journal: eLife
In common: statannotations, Keras, Numba, 12 other tools
[6] doi:10.1371/journal.pcbi.1014571 [code]
SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: Keras, Numba, TensorFlow, 11 other tools, 2 references
[7] doi:10.1038/s41467-026-72437-1 [code]
High-speed whole-brain imaging in Drosophila.
Journal: Nature communications
In common: OpenCV, scikit-image, h5py, 7 other tools, drosophila, systems, 4 references
[8] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Keras, UMAP, Numba, 12 other tools
[9] doi:10.1038/s41467-026-72152-x [code]
Centralized brain networks controlling antennal grooming coordination.
Journal: Nature communications
In common: TensorFlow, OpenCV, h5py, 7 other tools, drosophila, systems, 4 references
[10] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: xarray, Numba, OpenCV, 11 other tools, systems

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.