OSCR

Cheese3D enables sensitive detection and analysis of whole-face movement in mice.

Code ↔ Paper

19 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 19 matches
  1. [1] § Results › Uncovering underlying physiology from external facial movements ↔ paper/fig3-part1-cheese3d-general-anesthesia-eeg.ipynb, lines 553–643 · score 0.88 · EEG spectrogram, vertical scale bar, 0.2–1 Hz, EEG frequency band, FFT window, nose bulge volume
  2. [2] § Methods › Analysis of chewing kinematics ↔ paper/fig4-part1-chewing-whole-face-kinematics.ipynb, lines 681–738 · score 0.86 · find_peaks, median filtering, linearly interpolating, Kneedle, envelope, mouth area
  3. [3] § Methods › Analysis of in vivo electrical stimulation and electrophysiological recording ↔ fig5-part2-cheese3d-synchronized-electrophysiology.ipynb, lines 1110–1204 · score 0.82 · cyclic shuffling, corresponding lag, cross correlated, firing rates, facial movements, facial feature
  4. [4] § Methods › Anatomical-based interpretable feature selection ↔ packages/cheese3d/cheese3d/anatomy.py, lines 309–335 · score 0.81 · convex hull, right pad side, right pad top, nose bottom, axis, volume
  5. [5] § Results › Linking facial movement to motor control machineries using synchronized Cheese3D with electrophysiology ↔ fig5-part2-cheese3d-synchronized-electrophysiology.ipynb, lines 1110–1204 · score 0.81 · confidence interval, peak correlation, shuffled spike, cross correlation, firing rate, facial movements
  6. [6] § Methods › Analysis of in vivo electrical stimulation and electrophysiological recording ↔ fig5-part3-prediction-of-neural-activity-from-cheese3d.ipynb, lines 1113–1255 · score 0.80 · neural activity, cross validated, populationglm, intercept, GLMs, chunk
  7. [7] § Methods › Analysis of kinematics during anesthesia ↔ paper/fig3-part1-cheese3d-general-anesthesia-eeg.ipynb, lines 853–954 · score 0.79 · squared error, cross validation, eye height, nose bulge, PolynomialFeatures, ear angle
  8. [8] § Results › Reduction of tracking noise enables precise measurement of subtle and transient movements across facial regions ↔ paper/fig2-cheese3d-jitter-analysis.ipynb, lines 705–759 · score 0.75 · post triangulation, Jitter comparison, lateralized facial features, midline features, single mouse, camera view
  9. [9] § Results › Uncovering underlying physiology from external facial movements ↔ paper/fig3-part2-prediction-of-eeg-from-facial-features.ipynb, lines 1068–1157 · score 0.65 · sub delta, nose bulge volume, linear model, eye height, ear angle, VAR
  10. [10] § Results › Reduction of tracking noise enables precise measurement of subtle and transient movements across facial regions ↔ supfig5-cheese3d-jitter.ipynb, lines 925–1025 · score 0.63 · Whisker pad bulge, percentile jitter, Eye area, Nose bulge, Ear angle, facial region
  11. [11] § Results › Uncovering underlying physiology from external facial movements ↔ paper/fig3-part1-cheese3d-general-anesthesia-eeg.ipynb, lines 1155–1221 · score 0.61 · sub delta frequency, EEG power, frequency band power, theta, anesthetic, injection
  12. [12] § Results › Uncovering underlying physiology from external facial movements ↔ paper/fig3-part1-cheese3d-general-anesthesia-eeg.ipynb, lines 1223–1287 · score 0.56 · general anesthesia, movement raster, jitter threshold, facial movement, poses, vertical
  13. [13] § Methods › Video capture, synchronization and 3D calibration system ↔ packages/cheese3d/cheese3d/project.py, lines 478–557 · score 0.56 · square side length, board, ChArUco, pipeline, Anipose, triangulation
  14. [14] § Results › Uncovering underlying physiology from external facial movements ↔ paper/fig3-part2-prediction-of-eeg-from-facial-features.ipynb, lines 1377–1447 · score 0.54 · sub delta, EEG power, theta, latent, variance, predict
  15. [15] § Results › Linking facial movement to motor control machineries using synchronized Cheese3D with electrophysiology ↔ fig5-part3-prediction-of-neural-activity-from-cheese3d.ipynb, lines 1662–1774 · score 0.53 · explained variance, neural activity, ear angle, Poisson, predict, mice
  16. [16] § Results › Reduction of tracking noise enables precise measurement of subtle and transient movements across facial regions ↔ packages/cheese3d/cheese3d/interactive.py, lines 1–84 · score 0.53 · pose tracking, keypoint tracking, tool, Cheese3D, triangulation, training
  17. [17] § Methods › Analysis of kinematics during anesthesia ↔ paper/fig3-part2-prediction-of-eeg-from-facial-features.ipynb, lines 943–1027 · score 0.53 · fold cross validation, latent, smooth, predicted, EEG, linear
  18. [18] § Methods › Analysis of kinematics during anesthesia ↔ paper/fig3-part2-prediction-of-eeg-from-facial-features.ipynb, lines 1159–1243 · score 0.52 · latent model, eye height, nose bulge, ear angle, facial features, predict
  19. [19] § Methods › Analysis of chewing kinematics ↔ paper/fig4-part1-chewing-whole-face-kinematics.ipynb, lines 1178–1217 · score 0.52 · peak cross correlation, phases, mastication, ingestion, mouth area, chewing

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 · 1,578 lines · 65 KB · MIT · 4 matches

  1. # %% [markdown]
  2. # # FIGURE 3 (Feb 2025 submission)
  4. # %% [markdown]
  5. # This notebook gathers the code to make all the panels in Figure 3. In order to run this notebook, you need the following anipose projects:
  6. # 1. 20231013-long-anes-rig2
  7. # 2. 202408-eeg-emg-all
  9. # %% [markdown]
  10. # ## Prep the notebook
  12. # %% [markdown]
  13. # ### Load libraries
  15. # %%
  16. %load_ext autoreload
  17. %autoreload 2
  19. # Update path as if notebook was run from top-level repo directory
  20. import os
  21. import sys
  23. pwd = %pwd
  24. if pwd.endswith('notebooks'):
  25. sys.path.insert(0, os.path.abspath('..'))
  26. new_pwd = os.path.abspath(f"{pwd}/..")
  27. %cd {new_pwd}
  29. # %%
  30. import numpy as np
  31. import pandas as pd
  32. import matplotlib.pyplot as plt
  33. import matplotlib.ticker as ticker
  34. import seaborn as sns
  35. import itertools
  36. import pickle
  37. import json
  38. from glob import glob
  39. from numpy import newaxis as na
  40. from datetime import datetime
  41. from pathlib import Path
  42. from functools import reduce
  43. import scipy.stats as stats
  44. from scipy.stats import pearsonr, ttest_rel, wilcoxon
  45. from scipy.signal import convolve, spectrogram, sosfilt, iirdesign
  46. from sklearn.preprocessing import PolynomialFeatures
  47. from sklearn.model_selection import GridSearchCV
  48. from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
  49. from matplotlib.colors import ListedColormap
  50. from datetime import datetime
  52. from fepipeline.features.landmarks import read_3d_data
  53. from fepipeline.anatomy import compute_measurements_df
  54. from labutils.utils import maybe
  55. from labutils.plotting import (sns_setup,
  56. landmark_cmap,
  57. measurements_cmap,
  58. save_figure)
  60. # %% [markdown]
  61. # ### Global variables
  63. # %%
  64. sns_setup(palette="colorblind", font="sans-serif")
  65. OUTPUT_DIR = "./figures/Figure3-EEG"
  66. os.makedirs(OUTPUT_DIR, exist_ok=True)
  67. TODAY = datetime.today().date()
  69. # size in frames of moving average filter window for seeing slow drift in behavior
  70. FILTER_WINDOW = 1000
  72. # Start of anesthesia video after injection (in seconds)
  73. ANES_START_OFFSETS = {
  74. ("20240822", "B47"): 4 * 60 + 22,
  75. ("20240822", "B48"): 3 * 60 + 50,
  76. ("20240912", "B47"): 4 * 60 + 17,
  77. ("20240912", "B48"): 4 * 60 + 17,
  78. ("20240913", "B47"): 6 * 60 + 30, #anes video only
  79. ("20240913", "B48"): 3 * 60 + 33, #anes video only
  80. ("20240917", "B47"): 4 * 60 + 40,
  81. ("20240917", "B48"): 4 * 60 + 3,
  82. ("20240918", "B47"): 3 * 60 + 20,
  83. ("20240918", "B48"): 3 * 60 + 45,
  84. ("20240924", "B53"): 3 * 60 + 5, #anes video only
  85. ("20240925", "B53"): 5 * 60 + 56,
  86. ("20240926", "B53"): 4 * 60 + 3, #anes video only
  87. ("20240927", "B53"): 3 * 60 + 0, #anes video only
  88. }
  90. REGION_ORDER = ["eye(left)", "eye(right)",
  91. "ear(left)", "ear(right)",
  92. "nose", "cheek", "mouth"]
  94. DATA_CACHE = "measurements-data-cache-2024"
  95. os.makedirs(DATA_CACHE, exist_ok=True)
  96. MEAS_DATA_CACHE = os.sep.join([DATA_CACHE, "eeg-slow-drift.pkl"])
  97. ALL_DATA_CACHE = os.sep.join([DATA_CACHE, "eeg-slow-drift-eeg-meas.pkl"])
  98. MODEL_RESULTS_CACHE = os.sep.join([DATA_CACHE, "long-anes-slow-drift-results.pkl"])
  100. # Data cache for time-prediction models using Lasso - Alpha optimization
  101. # Using facial features
  102. LASSO_MODEL_FF = os.sep.join([DATA_CACHE, f"{TODAY}-lasso-model-ff.pkl"])
  103. # Using eeg features
  104. LASSO_MODEL_EEG = os.sep.join([DATA_CACHE, f"{TODAY}-lasso-model-eeg.pkl"])
  106. ANIPOSE_BASE = 'anipose-projects/202408-eeg-emg-all'
  107. COORDINATE_PATHS = {}
  108. key_cols = ('date', 'mouse', 'condition')
  109. for p in Path(ANIPOSE_BASE).glob('*/pose-3d/*.csv'):
  110. date, mouse, *_ = p.name.split('_')
  111. if date in ["20240816", "20240822"]:
  112. continue
  113. condition = 'awake' if 'awake' in p.name else 'anes'
  114. COORDINATE_PATHS[(date, mouse, condition)] = p
  115. data_keys = list(COORDINATE_PATHS.keys())
  116. anes_data_keys = [d for d in data_keys if d[2] != 'awake']
  118. # %%
  119. def compute_sample_rate(timestamps):
  120. start_time = datetime.strptime(timestamps[0], "%H:%M:%S:%f")
  121. times = np.array([(datetime.strptime(t, "%H:%M:%S:%f") - start_time).total_seconds()
  122. for t in timestamps])
  124. return 1 / np.mean(np.diff(times))
  126. def read_signal(filename, sample_rate):
  127. df = pd.read_csv(filename, sep="\t", names=["timestamp", "signal"])
  128. _sample_rate = compute_sample_rate(df["timestamp"].values)
  129. if abs(sample_rate - _sample_rate) / sample_rate > 0.1:
  130. return df["signal"].values, _sample_rate
  131. else:
  132. return df["signal"].values, sample_rate
  134. # Remove NaNs at the beginning of the EEG recording and define new recording start
  135. EEG_STARTS = os.sep.join([DATA_CACHE, 'eeg-anes-start.pkl'])
  136. if os.path.exists(EEG_STARTS):
  137. with open(EEG_STARTS, "rb") as fio:
  138. ANES_START_OFFSETS_EEG = pickle.load(fio)
  139. else:
  140. ANES_START_OFFSETS_EEG = {}
  141. for date, mouse, cond in data_keys:
  142. folder = glob(f"ephys-data/{date}_*{mouse}*_EEG-EMG-rec_rig2")
  143. with open(os.sep.join([folder[0], f"{date}_{mouse}_{cond}.align.json"])) as f:
  144. alignment = json.load(f)
  145. lag_time = alignment["lag_time"]
  146. sample_rate = alignment["sample_rate"]
  147. for name, units in (("eeg", "V"), ("emg", "V"), ("temp", "C")):
  148. filename = os.sep.join([folder[0], f"{date}_{mouse}_{cond}_{name}.txt"])
  149. signal, _sample_rate = read_signal(filename, sample_rate)
  150. # if name in ['emg', 'eeg']:
  151. signal[np.isnan(signal)] = 0 #added by inm - assign the nan to 0s then save signal from first non-zero
  152. signal_start = signal.nonzero()[0][0]
  153. signal = signal[signal_start:]
  154. if name == 'eeg' and signal_start and cond == 'anes':
  155. ANES_START_OFFSETS_EEG[date, mouse] = ANES_START_OFFSETS[date, mouse] + int(signal_start/_sample_rate)
  157. # %% [markdown]
  158. # ### Define colormaps
  160. # %%
  161. _, LANDMARK_CMAP = landmark_cmap()
  162. _, MEASUREMENT_CMAP = measurements_cmap()
  164. CONTROL_CMAP = sns.color_palette([sns.color_palette("colorblind")[2],
  165. sns.color_palette("colorblind")[-3],
  166. sns.color_palette("colorblind")[3]])
  168. FREQ_CMAP = "jet"
  170. # %%
  171. EEG_CONTROL_CMAP = sns.color_palette([sns.color_palette("colorblind")[0],
  172. sns.color_palette("colorblind")[7],
  173. sns.color_palette("colorblind")[0],
  174. sns.color_palette("colorblind")[0]])
  175. EEG_CONTROL_CMAP
  177. # %%
  178. CONTROL_CMAP_V2 = sns.color_palette([sns.color_palette("colorblind")[3],
  179. sns.color_palette("colorblind")[-3],
  180. sns.color_palette("colorblind")[2],
  181. sns.color_palette("colorblind")[0],
  182. sns.color_palette("colorblind")[1],
  183. sns.color_palette("colorblind")[-1],
  184. sns.color_palette("colorblind")[4],
  185. sns.color_palette("colorblind")[-2]])
  186. CONTROL_CMAP_V2
  188. # %% [markdown]
  189. # ### Define functions
  191. # %%
  192. def z_score(raw_data):
  193. mean = np.mean(raw_data, 0, keepdims=True)
  194. sigma = np.std(raw_data, 0, keepdims=True) #1
  195. sigma[sigma == 0] = 1
  196. return (raw_data-mean)/sigma
  198. def un_zscore(raw_data, z_score_data): #Undo z-scoring for plotting as a sanity check - there is no deformation
  199. mean = np.mean(raw_data, 0, keepdims=True) #0
  200. sigma = np.std(raw_data, 0, keepdims=True) #1
  201. sigma[sigma == 0] = 1
  202. return z_score_data*sigma + mean
  204. # %%
  205. def moving_avg_filter(signal, window = 3, stride = 1, truncate = True):
  206. kernel = np.ones(window) / window
  207. result = convolve(signal, kernel, mode="same")[window:-window]
  209. return result[::stride]
  211. # %%
  212. test = moving_avg_filter(np.arange(12), window = 1, stride = 2)
  213. print(test)
  215. # %%
  216. def plot_spectrum_helper(dates, mice, conditions, spectrums, lag_times, color,
  217. min_power = None, max_power = None):
  218. single = (len(spectrums.values) == 1)
  219. if single:
  220. if min_power is None:
  221. min_power = np.min(spectrums.values[0][-1])
  222. if max_power is None:
  223. max_power = np.max(spectrums.values[0][-1])
  224. else:
  225. if min_power is None:
  226. min_power = min(*(np.min(s) for _, _, s in spectrums.values))
  227. if max_power is None:
  228. max_power = max(*(np.max(s) for _, _, s in spectrums.values))
  229. awake = (conditions == "awake")
  230. # print(awake)
  231. if sum(awake):
  232. fs, ts, power = spectrums[awake].values[0]
  233. awake_end = ts[-1]
  234. else:
  235. awake_end = 0
  236. fs, ts, power = spectrums[~awake].values[0]
  237. anes_start = ANES_START_OFFSETS[(dates.values[0], mice.values[0])]
  238. ts = ts + lag_times[~awake].values[0] + anes_start #+ awake_end
  239. # print(power)
  240. p = plt.pcolormesh(ts / 60, fs, power,
  241. vmin=min_power, vmax=max_power, cmap=FREQ_CMAP)
  244. # %%
  245. def compute_power_bands(fs, power):
  246. sub_delta_idx = np.where(fs < 1)[0]
  247. sub_delta_power = np.mean(power[sub_delta_idx], axis=0)
  249. delta_idx = np.where((fs >= 1) & (fs <= 4))[0]
  250. delta_power = np.mean(power[delta_idx], axis=0)
  252. theta_idx = np.where((fs >= 5) & (fs <= 10))[0]
  253. theta_power = np.mean(power[theta_idx], axis=0)
  255. # slow_delta_idx = np.where(fs <= 4)[0]
  256. # slow_delta_power = np.mean(power[slow_delta_idx], axis=0)
  258. return sub_delta_power, delta_power, theta_power #, slow_delta_power
  260. # %%
  261. def generate_lagged_data(data, lags = 1, bias = True, initial_bias = True):
  262. assert data.ndim == 2 or data.ndim == 1
  263. # make the data a matrix
  264. if data.ndim == 1:
  265. data = np.expand_dims(data, axis=1)
  266. # make shifted copies of data
  267. lagged = np.lib.stride_tricks.sliding_window_view(data, lags, axis=0)
  268. lagged = np.reshape(lagged, (data.shape[0] - lags + 1, -1))
  269. # pad in the initial lags
  270. if lags > 1:
  271. pad = np.stack([np.pad(data[:(lags - i)], ((i, 0), (0, 0))).T
  272. for i in range(lags - 1, 0, -1)])
  273. pad = np.reshape(pad, (lags - 1, -1))
  274. lagged = np.concatenate([pad, lagged], axis=0)
  275. # always have initial time point
  276. if initial_bias:
  277. lagged = np.concatenate([np.tile(data[0], (lagged.shape[0], 1)), lagged], axis=1)
  278. # add bias term
  279. if bias:
  280. lagged = np.pad(lagged, ((0, 0), (1, 0)), constant_values=1)
  282. return lagged
  284. # Uncomment prints for examples:
  285. x = np.arange(10)
  286. y = np.stack([x, x + 10], axis=1)
  287. # print(generate_lagged_data(x, 5))
  288. # print(generate_lagged_data(y, 3))
  290. # %%
  291. def generate_poly_feats(data, degree = 2, interaction = True):
  292. assert data.ndim == 2 or data.ndim == 1
  293. # make the data a matrix
  294. if data.ndim == 1:
  295. data = np.expand_dims(data, axis=1)
  296. if interaction == True:
  297. return PolynomialFeatures(degree).fit_transform(data)
  298. else:
  299. polyfeats = np.vstack((np.ones(len(data)), data.T)).T
  300. for n in (np.arange(1,degree)+1):
  301. polyfeats = np.hstack((polyfeats, data**n))
  303. return polyfeats
  305. # Uncomment prints for examples:
  306. z = np.stack([x, x + 10, x + 1], axis=1)
  307. print(generate_lagged_data(x, bias=False))
  308. print(generate_poly_feats(generate_lagged_data(z, bias=False)))
  309. print(generate_poly_feats(generate_lagged_data(x, bias=False), 3))
  311. # %%
  312. def get_features(measures, degree = 2, lags = 1, interaction = True, initial_bias = True):
  313. features = generate_lagged_data(measures, lags, bias=False, initial_bias = initial_bias)
  314. features = generate_poly_feats(features, degree, interaction)
  316. return features
  318. # %%
  319. def plot_predictions(ts, xs, color, label, FILTER_WINDOW = 1000, dash_line = True, ax = None):
  320. ts = ts.values[0][FILTER_WINDOW:-FILTER_WINDOW]
  321. xs = moving_avg_filter(xs.values[0], FILTER_WINDOW)
  322. if ax is None:
  323. sns.lineplot(x=ts, y=xs, color=color, label=label)
  324. if dash_line:
  325. sns.lineplot(x=ts, y=ts, color="black", linestyle="dashed")
  326. else:
  327. sns.lineplot(x=ts, y=xs, color=color, label=label, ax = ax)
  328. if dash_line:
  329. sns.lineplot(x=ts, y=ts, color="black", linestyle="dashed", ax = ax)
  331. # %%
  332. def format_feature_ticks(feature):
  333. text = feature.get_text()
  334. if text == "whole-face":
  335. return "whole\nface"
  336. if ' ' in text:
  337. parts = text.split(" ")
  338. return'\n'.join(parts)
  339. else:
  340. parts = text.split(", ")
  341. return ',\n'.join(parts)
  343. # %%
  344. def test_mean_per_run(results, n_runs, feat):
  345. mean_test = []
  346. for r in np.arange(n_runs):
  347. mean_test.append([r, feat, np.mean(results.query("set == 'Test' & run == @r")["RMSE"].values)])
  349. return pd.DataFrame(mean_test, columns = ["run",
  350. "features",
  351. "mean_RMSE"])
  353. # %% [markdown]
  354. # ### Load data
  356. # %%
  357. if os.path.exists(MEAS_DATA_CACHE):
  358. meas_df = pd.read_pickle(MEAS_DATA_CACHE)
  359. else:
  360. coord_data = {k: read_3d_data(v.parent.parent.as_posix())
  361. for k, v in COORDINATE_PATHS.items()}
  362. meas_df = compute_measurements_df(coord_data, key_columns=key_cols)
  363. meas_df = meas_df.assign(sample_rate=100.0, lag_time=0.0)
  364. pd.to_pickle(meas_df, MEAS_DATA_CACHE)
  365. meas_df = meas_df.drop(columns=["measurement_value", "std", "count"])
  366. meas_df
  368. # %%
  369. if os.path.exists(ALL_DATA_CACHE):
  370. data = pd.read_pickle(ALL_DATA_CACHE)
  371. else:
  372. eeg_rows = []
  373. for date, mouse, cond in data_keys:
  374. folder = glob(f"ephys-data/{date}_*{mouse}*_EEG-EMG-rec_rig2")
  375. if len(folder) != 1:
  376. print(f"Found too many (or no) sources={folder} for {date=}, {mouse=}, {cond=}")
  377. continue
  378. with open(os.sep.join([folder[0], f"{date}_{mouse}_{cond}.align.json"])) as f:
  379. alignment = json.load(f)
  380. lag_time = alignment["lag_time"]
  381. sample_rate = alignment["sample_rate"]
  382. for name, units in (("eeg", "V"), ("emg", "V"), ("temp", "C")):
  383. filename = os.sep.join([folder[0], f"{date}_{mouse}_{cond}_{name}.txt"])
  384. signal, _sample_rate = read_signal(filename, sample_rate)
  385. # if name in ['emg', 'eeg']:
  386. signal[np.isnan(signal)] = 0 #added by inm - assign the nan to 0s then save signal from first non-zero
  387. signal_start = signal.nonzero()[0][0]
  388. signal = signal[signal_start:]
  389. if name == 'eeg' and signal_start:
  390. ANES_START_OFFSETS[date, mouse] = ANES_START_OFFSETS[date, mouse] + signal_start/_sample_rate
  391. eeg_rows.append([date, mouse, cond,
  392. "ephys", name, units, "ephys",
  393. signal, _sample_rate, lag_time])
  395. data = pd.concat([meas_df, pd.DataFrame(eeg_rows, columns=meas_df.columns)])
  396. pd.to_pickle(data, ALL_DATA_CACHE)
  398. data
  400. # %%
  401. from copy import deepcopy
  402. def compute_spectrogram(row):
  403. signal = row["timeseries"]
  404. # if np.sum(np.isnan(signal))>0:
  405. # print("signal", np.sum(np.isnan(signal)))
  406. # signal[np.isnan(signal)] = 0 #added by inm - assign the nan to 0s so that they do not propagate when computing spectrogram - to be improved
  408. eeg_filt = iirdesign(50, 55, 1, 50, fs=row["sample_rate"], output="sos")
  409. fs, ts, power = spectrogram(sosfilt(eeg_filt, signal),
  410. fs=row["sample_rate"],
  411. nperseg= 5000,#2000, #originally 2048 ; possibly change to 2050 or 2000 (divisible by 10) - number of samples in each fft segment
  412. # noverlap=512,
  413. scaling="spectrum",
  414. mode="magnitude")
  415. # if np.sum(np.isnan(power)):
  416. # print("power", np.sum(np.isnan(power))) #nans are not here but appear in the plot
  417. row = row.copy()
  418. row["measurement_name"] = "eeg-power"
  419. row["measurement_units"] = "V^2"
  420. row["timeseries"] = (fs, ts, power)
  422. return row
  424. eeg_df = deepcopy(data.query("measurement_name == 'eeg'")) #brand new copy so that I can pass by object (instead of by reference) and do not alter the original data
  425. data_v2 = pd.concat([data, eeg_df.apply(compute_spectrogram, axis=1)])
  426. data_v2
  428. # %%
  429. data_and_bands = deepcopy(data_v2)
  430. eeg_window = 8 #20
  431. band_names = ['subdelta', 'delta', 'theta']
  432. for d, m, c in anes_data_keys:
  433. fs, ts, power = data_and_bands.query("mouse == @m & date == @d & condition == @c & measurement_name == 'eeg-power'")["timeseries"].values[0]
  434. bands = compute_power_bands(fs, power)
  435. for i_b, _b in enumerate(bands):
  436. bands_dict = {'date': d,
  437. 'mouse': m,
  438. 'condition': c,
  439. 'measurement_group': 'ephys',
  440. 'measurement_name': band_names[i_b],
  441. 'measurement_units': 'V^2',
  442. 'measurement_type': 'ephys',
  443. #Save power band values and actual time together
  444. 'timeseries': (z_score(moving_avg_filter(_b, window = eeg_window, stride = 1)),
  445. ts[eeg_window:-eeg_window]),
  446. 'sample_rate': data_and_bands.query("mouse == @m & date == @d & condition == @c & measurement_name == 'eeg-power'")["sample_rate"].values[0],
  447. 'lag_time': data_and_bands.query("mouse == @m & date == @d & condition == @c & measurement_name == 'eeg-power'")["sample_rate"].values[0]}
  448. data_and_bands = pd.concat([data_and_bands, pd.DataFrame.from_records([bands_dict])])
  450. data_and_bands
  453. # %%
  454. print('Example mouse B48, date 20240917, condition anes.')
  455. print('Number of raw EEG samples: ', data.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg' & condition == 'anes'")["timeseries"].values[0].shape)
  456. print('Sample rate: ', data.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg' & condition == 'anes'")["sample_rate"].values[0])
  457. print('Number of power EEG samples: ', data_v2.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg-power' & condition == 'anes'")["timeseries"].values[0][1].shape)
  458. print('The EEG signal was downsampled by: ',
  459. data.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg' & condition == 'anes'")["timeseries"].values[0].shape[0]/
  460. data_v2.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg-power' & condition == 'anes'")["timeseries"].values[0][1].shape[0])
  461. new_sr = data.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg' & condition == 'anes'")["sample_rate"].values[0]/(data.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg' & condition == 'anes'")["timeseries"].values[0].shape[0]/
  462. data_v2.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg-power' & condition == 'anes'")["timeseries"].values[0][1].shape[0])
  463. print('Therefore the new sampling rate is: ', new_sr)
  464. print('And first three timepoints of EEG power are found at: ', data_v2.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg-power' & condition == 'anes'")["timeseries"].values[0][1][0:3], 'seconds')
  466. print('\nDownsampling by the spectrum function is more or less constant across mice and sessions:')
  467. print(data.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg' & condition == 'anes'")["timeseries"].values[0].shape[0]/
  468. data_v2.query("mouse == 'B48' & date == '20240917' & measurement_name == 'eeg-power' & condition == 'anes'")["timeseries"].values[0][1].shape[0])
  469. print(data.query("mouse == 'B47' & date == '20240918' & measurement_name == 'eeg' & condition == 'anes'")["timeseries"].values[0].shape[0]/
  470. data_v2.query("mouse == 'B47' & date == '20240918' & measurement_name == 'eeg-power' & condition == 'anes'")["timeseries"].values[0][1].shape[0])
  472. # %%
  473. cam_sr = 100
  474. eeg_power_sr = {}
  475. ds_factor = {}
  476. for date, mouse, _ in anes_data_keys:
  477. eeg_power_sr[mouse, date] = data_v2.query("mouse == @mouse & date == @date & measurement_name == 'eeg' & condition == 'anes'")["sample_rate"].values[0]/(
  478. data_v2.query("mouse == @mouse & date == @date & measurement_name == 'eeg' & condition == 'anes'")["timeseries"].values[0].shape[0]/
  479. data_v2.query("mouse == @mouse & date == @date & measurement_name == 'eeg-power' & condition == 'anes'")["timeseries"].values[0][1].shape[0])
  480. ds_factor[mouse, date] = int(np.floor(cam_sr / eeg_power_sr[mouse, date])) # camera sample rate / eeg power sample rate
  482. ds_factor
  484. # %% [markdown]
  485. # ## Figure 3c
  487. # %% [markdown]
  488. # (notebook: 2024-12-05-long-anes-three-way-predictions-eeg.ipynb)
  490. # %% [markdown]
  491. # Overlaid z-scored ear angle, eye height, and nose bulge volume traces from 12 sessions (n = 3 mice, four sessions per mouse. Vertical scale bar indicates one standard deviation).
  493. # %%
  494. # f_window = 1000 #1000
  495. features = ['ear-angle-left', 'eye-height-left', 'nose-bulge-volume']
  496. fig, axs = plt.subplots(nrows=len(features),
  497. ncols=1,
  498. figsize=(3*3, 4*3),
  499. dpi=200,
  500. constrained_layout=False)
  501. meas_summary = {}
  502. time_summary = {}
  503. for row, l in enumerate(features):
  504. meas_summary[l] = []
  505. time_summary[l] = []
  506. y_labels = []
  508. for sess in anes_data_keys:
  509. mouse = sess[1]
  510. date = sess[0]
  511. if sess in anes_data_keys:
  512. anes_start = np.max(list(ANES_START_OFFSETS.values())) - ANES_START_OFFSETS[date, mouse] #in seconds - Take the same start for all of them
  513. ys = data_v2.query(
  514. "mouse == @mouse & "
  515. "date == @date & "
  516. "condition == 'anes' & "
  517. "measurement_name == @l"
  518. )["timeseries"].values[0]
  519. anes_end = int(np.floor((45 * 60 - ANES_START_OFFSETS[date, mouse]))) * 100 # Trim at 45 minutes since anesthesia injection
  520. ys = ys[int(np.floor(anes_start)*100):anes_end]
  522. time = (np.arange(len(ys))/100 + np.max(list(ANES_START_OFFSETS.values())))/ 60 #Note all will start at 5 min and camera sample rate is 100
  523. time = time[FILTER_WINDOW:-FILTER_WINDOW]
  524. # FILTER AND Z SCORED DATA PER MEASUREMENT
  525. ys = moving_avg_filter(ys, FILTER_WINDOW)
  526. ys = z_score(ys)
  527. meas_summary[l].append(ys)
  528. time_summary[l].append(time)
  529. sns.lineplot(x = time, y = ys,
  530. color = CONTROL_CMAP_V2[0], linewidth = 0.5, alpha = 0.5,
  531. ax = axs[row])
  532. sns.lineplot(x = time, y = np.mean(meas_summary[l], axis = 0),
  533. color = CONTROL_CMAP_V2[0], linewidth = 2.5,
  534. ax = axs[row])
  536. axs[row].sharex(axs[0])
  537. axs[row].set_ylabel(('\n').join(l.split('-')), rotation = 0)
  538. axs[row].yaxis.set_label_coords(-0.1,0.5)
  539. if row != (len(features)-1):
  540. axs[row].tick_params("x", bottom=False, labelbottom=False)
  541. sns.despine(ax=axs[row], bottom=True)
  544. axs[row].set_xlabel('Time (s)')
  545. axs[row].set_xticks(np.arange(0, 45, 10))
  546. sns.despine(ax=axs[row])
  548. fig.tight_layout()
  550. fname = "3c-facial-feat-3subset-allmice-10s-avgfilter"
  551. fig.savefig(os.sep.join([OUTPUT_DIR, f"{fname}.svg"]), bbox_inches="tight")
  553. # %% [markdown]
  554. # ## Figure 3d
  556. # %% [markdown]
  557. # (notebook: 2024-08-30-eeg-emg-long-anes.ipynb)
  559. # %% [markdown]
  560. # Simultaneously recorded Cheese3D features (top, showing moving average over a 10 sec window; vertical scale bars: ear: 2◦, eye: 0.1mm, nose: 0.5mm3), EEG spectrogram (middle, 5sec FFT window), and power of EEG frequency bands (bottom, showing subdelta: 0.2 Hz to 1 Hz, delta: 1 Hz to 4 Hz, and theta: 5 Hz to 10 Hz bands; vertical scale bars indicate one standard deviation) for an example session.
  562. # %%
  563. max_power = np.max(data_v2.query("measurement_name == 'eeg-power'").apply(
  564. lambda x: np.max(x["timeseries"][-1]),
  565. axis=1
  566. ).values)
  567. min_power = np.min(data_v2.query("measurement_name == 'eeg-power'").apply(
  568. lambda x: np.min(x["timeseries"][-1]),
  569. axis=1
  570. ).values)
  572. def plot_aligned_helper(dates, mice, conditions, style, timeseries, lags, sample_rate, color):
  573. anes_start = ANES_START_OFFSETS[(dates.values[0], mice.values[0])] #in seconds
  574. anes_end = -1 #np.floor((45 * 60 - ANES_START_OFFSETS[(dates.values[0], mice.values[0])])) # in seconds - Trim at 45 minutes since anesthesia injection
  575. awake = (conditions == "awake")
  576. if style.values[0] == "eeg-power":
  577. plot_spectrum_helper(dates, mice, conditions, timeseries, lags, color,
  578. min_power=min_power, max_power=max_power)
  579. elif style.values[0] in band_names:
  580. _b = timeseries[~awake].values[0][0] #this has already been filtered and zscored (see cell above)
  581. ts = timeseries[~awake].values[0][1]
  582. sns.lineplot(x=(ts[:anes_end]+anes_start)/60,
  583. y=_b[:anes_end])
  584. else:
  585. ys = timeseries[~awake].values[0]
  586. ys = ys[:anes_end]
  587. xs = np.arange(len(ys)) / sample_rate[~awake].values[0]
  588. xs = xs + anes_start #+ awake_end
  589. sns.lineplot(x=xs[FILTER_WINDOW:-FILTER_WINDOW] / 60,
  590. y=moving_avg_filter(ys, FILTER_WINDOW),
  591. color=CONTROL_CMAP[2],
  592. label="anes")
  594. measurements = ["ear-angle-left",
  595. "eye-height-left",
  596. "nose-bulge-volume",
  597. "eeg-power",
  598. "theta",
  599. "delta",
  600. "subdelta",
  601. ]
  602. for d in set(data_and_bands["date"]):
  603. if d == '20240822':
  604. continue
  605. g = sns.FacetGrid(data_and_bands.query("date == @d & "
  606. "measurement_name in @measurements"),
  607. row="measurement_name", col="mouse",
  608. row_order=measurements,
  609. sharey=False,
  610. sharex=False,
  611. # xlim=(0, 70),
  612. aspect=2,
  613. height=2)
  614. g.map(plot_aligned_helper,
  615. "date", "mouse", "condition",
  616. "measurement_name", "timeseries", "lag_time", "sample_rate")
  617. cbar_ax = g.figure.add_axes([1.015, 0.44, 0.015, 0.12]) #[1.015, 0.19, 0.015, 0.12]
  618. plt.colorbar(cax=cbar_ax)
  619. g.figure.suptitle(f"Anesthetized Spectrogram ({d})")
  620. g.set_xlabels("Time (min)")
  621. g.set_ylabels("")
  622. sns.despine(g.figure)
  623. for ax in g.axes.flat:
  624. # if ax.get_xbound()[1] > 65:
  625. # ax.set_xticks(np.arange(0, 71, 10))
  626. # else:
  627. ax.set_xticks(np.arange(0, 45, 10))
  628. ax.set_xlim(0, 45)
  629. for ax in g.axes[-4, :]: # Limit for the EEG spectrogram
  630. ax.set_ylim(0, 20)
  631. for ax in g.axes[:-1, :].flat: # Remove x-ticks except in bottom plot
  632. ax.set_xticks([])
  633. sns.despine(ax=ax, bottom=True)
  634. # for ax in g.axes[-1, :]:
  635. # axmin = min(ax.get_xlim()[0], ax.get_ylim()[0])
  636. # axmax = max(ax.get_xlim()[1], ax.get_ylim()[1])
  637. # ax.plot(np.linspace(axmin, axmax), np.linspace(axmin, axmax), color='black', linestyle="--")
  638. # ax.set_aspect('equal', anchor = 'SW', adjustable = 'box')
  640. g.set_titles(template="{col_name} ({row_name})")
  641. g.tight_layout()
  642. fname = f"3d-fe-and-eeg-exemplar-{d}"
  643. # g.figure.savefig(os.sep.join([OUTPUT_DIR, f"{fname}.svg"]), bbox_inches="tight")
  645. # %% [markdown]
  646. # ## Figure 3e
  648. # %% [markdown]
  649. # (notebook: 2024-12-05-long-anes-three-way-predictions-eeg)
  651. # %% [markdown]
  652. # Output from a quadratic model fit across mice predicting time since injection using the initial and current Cheese3D (orange) or EEG (blue) feature values relative to the dotted identity line.
  654. # %%
  655. EEG_MODEL_NAME = os.sep.join([DATA_CACHE, "lasso-model-mouse-sess-eeg-220runs-interact-stride-init-bias-45minend-Jan30-2ndfix.pkl"])
  656. FF_MODEL_NAME = os.sep.join([DATA_CACHE, "lasso-model-mouse-sess-220runs-interact-stride-init-bias-45minend-Jan30-2ndfix.pkl"])
  658. results_eeg_df_name = 'measurements-data-cache-2024/prediction-eeg-results-220runs-interact-stride-init-bias-45minend-Jan30-2ndfix.pkl'
  659. results_df_name = 'measurements-data-cache-2024/prediction-results-220runs-interact-stride-init-bias-45minend-Jan30-2ndfix.pkl'
  662. # %% [markdown]
  663. # Note that the names above correspond to the files used for the manuscript results. If you would like to use a different model, you will need to change the specified names.
  664. #
  665. # Comment the cell below if you would like to run the model code.
  667. # %%
  668. # These are the models and predictions used in the manuscript
  669. # EEG models and results
  670. if os.path.exists(EEG_MODEL_NAME):
  671. with open(EEG_MODEL_NAME, "rb") as fio:
  672. models_eeg = pickle.load(fio) # Model results
  673. if os.path.exists(results_eeg_df_name):
  674. results_eeg_df = pd.read_pickle(results_eeg_df_name) # Prediction results
  676. # Facial feature models
  677. if os.path.exists(FF_MODEL_NAME):
  678. with open(FF_MODEL_NAME, "rb") as fio:
  679. models_ff = pickle.load(fio)
  680. if os.path.exists(results_eeg_df_name):
  681. results_df = pd.read_pickle(results_df_name)
  683. # %% [markdown]
  684. # ### EEG model
  686. # %%
  687. # def run_regression(ts_df, anes_data_keys, features, train_idx, alpha_range, models = None):
  688. def run_regression(ts_df, anes_data_keys, alpha_range):
  690. # specify feature parameters
  691. # 20 filter window in samples ~ 40 seconds if nperseg = 2000;
  692. # 8 filter window in samples ~ 40 seconds if nperseg = 5000
  693. # refer to prints about eeg-power sample rate in cell above)
  694. eeg_window = 8
  696. # Polynomial features
  697. degree = 2 #degree of the polynomial
  698. lag = 1
  700. # define empty lists to generate data for model
  701. measures = []
  702. target = []
  704. indeces = []
  705. counter = 0
  706. for s_idx, dk in enumerate(anes_data_keys):
  707. d = dk[0]
  708. mouse = dk[1]
  710. # Select input features: eeg-power bands
  711. eegpower = ts_df.query("mouse == @mouse & date == @d & condition == 'anes' & measurement_name == 'eeg-power'")["timeseries"]
  712. fs, ts, power = eegpower.values[0]
  713. _bands = compute_power_bands(fs, power) #type(_bands) = tuple
  715. # Define actual time
  716. anes_start = ANES_START_OFFSETS_EEG[(d, mouse)]
  717. anes_end = int(np.floor((45 * 60 - anes_start) * eeg_power_sr[mouse, d])) # Trim at 45 minutes since anesthesia injection
  718. _time = ts[:anes_end]
  719. target.append((_time[eeg_window:-eeg_window] + anes_start) / 60) #lag is only for plotting comparisons against facial features?
  721. # Filter noise
  722. bands =[]
  723. for _b in _bands:
  724. bands.append(moving_avg_filter(_b[:anes_end], window = eeg_window, stride = 1))
  726. # Z-score and get features
  727. _measures = np.concatenate((z_score(bands[0][:,na]), #subdelta
  728. z_score(bands[1][:,na]), #delta
  729. z_score(bands[2][:,na])), 1) #theta
  730. _measures = get_features(_measures, degree, lag, interaction = True, initial_bias = True)
  731. _measures[:,0] = 1
  733. # Save indeces in a list of arrays corresponding to each session
  734. indeces.append(np.arange(_measures.shape[0]) + counter)
  735. counter += _measures.shape[0]
  736. measures.append(_measures)
  738. measures = np.concatenate(measures, axis=0)
  739. time = np.concatenate(target, axis=0)
  740. cv = []
  741. # Annotate indeces to separate into validation and train sets
  742. for k in range(len(indeces)):
  743. train_indeces = np.concatenate([indeces[i] for i in range(len(indeces)) if i != k])
  744. test_indeces = indeces[k] #Leave one session out for validation (a.k.a. test within the train set)
  745. cv.append((train_indeces, test_indeces))
  747. print(f"Running regression for (subdelta, delta, supdelta)...")
  749. # fit Lasso model and explore alpha space
  750. mdl = Lasso(fit_intercept=True)
  751. param_grid = {'alpha': alpha_range} #key must be the same name that is used in Lasso documentation
  752. model = GridSearchCV(mdl, param_grid, cv = cv, verbose = 0, n_jobs = 8, scoring = "neg_root_mean_squared_error") # Explore hyper-parameter space (only for alpha in this case)
  753. model.fit(measures[:,1:], time)
  755. return model
  758. # Explore alpha space below:
  759. alpha_range = np.logspace(-1, 2, 100)
  761. n_runs = 220 #This number should match the len(all_combos) below
  762. n_test = 3 #number of test sessions
  763. test_heldout =[]
  764. all_combos = list(itertools.combinations(range(len(anes_data_keys)), n_test)) #all possible combinations
  766. if os.path.exists(EEG_MODEL_NAME):
  767. print("Manuscript model has been loaded.")
  768. elif os.path.exists(LASSO_MODEL_EEG):
  769. print(f"{LASSO_MODEL_EEG} has been loaded.")
  770. with open(LASSO_MODEL_EEG, "rb") as fio:
  771. models_eeg = pickle.load(fio) # Model results
  772. else:
  773. models_eeg = []
  774. for r in np.arange(n_runs):
  775. test_heldout_idx = all_combos[r]
  776. _test_heldout = [d for idx, d in enumerate(anes_data_keys) if idx in test_heldout_idx]
  777. test_heldout.append(_test_heldout)
  778. anes_data_keys_train = [d for d in anes_data_keys if d not in _test_heldout]
  779. models_eeg.append(run_regression(data_v2, anes_data_keys_train, alpha_range))
  780. with open(LASSO_MODEL_EEG, "wb") as fio:
  781. pickle.dump(models_eeg, fio)
  785. # %%
  786. def get_results_df(data, model, anes_data_keys, test_heldout, n_run = 0):
  787. eeg_window = 8
  788. results = []
  790. for d in anes_data_keys:
  791. date = d[0]
  792. mouse = d[1]
  793. anes_start = ANES_START_OFFSETS_EEG[(date, mouse)] # in seconds
  794. anes_end = int(np.floor((45 * 60 - anes_start) * eeg_power_sr[mouse, date])) # Trim at 45 minutes since anesthesia injection
  796. # Raw eeg/power data
  797. eegpower = data.query("mouse == @mouse & date == @date & condition == 'anes' & measurement_name == 'eeg-power'")["timeseries"]
  798. fs, ts, power = eegpower.values[0]
  799. _bands = compute_power_bands(fs, power) #type(_bands) = tuple
  801. # Filter eeg power data
  802. bands =[]
  803. for _b in _bands:
  804. bands.append(moving_avg_filter(_b[:anes_end], window = eeg_window, stride = 1))
  806. # Z-score and get features
  807. measures = np.concatenate((z_score(bands[0][:,na]), #subdelta
  808. z_score(bands[1][:,na]), #delta
  809. z_score(bands[2][:,na])), 1) #theta
  810. test_ts = get_features(measures, degree = 2, lags = 1, interaction = True, initial_bias = True)
  811. test_ts[:,0] = 1
  813. # Define actual time
  814. time = ts[:anes_end]
  815. time = (time[eeg_window:-eeg_window] + anes_start) / 60 #lag is only for plotting comparisons against facial features?
  817. # predict time
  818. times_hat = model.predict(test_ts[:,1:])
  820. # Compute RMSE
  821. rmse = np.sqrt(np.mean((times_hat - time) ** 2))
  823. # Save results
  824. set_group = "Test" if (date, mouse, 'anes') in test_heldout else "Train"
  825. results.append([n_run, mouse, date, set_group, "(subdelta, delta, theta)", times_hat, time, rmse])
  827. results_df = pd.DataFrame(results, columns = ["run",
  828. "mouse",
  829. "date",
  830. "set",
  831. "features",
  832. "predicted time",
  833. "actual time",
  834. "RMSE"])
  835. return results_df
  837. _results_df = []
  838. PREDICTION_RESULTS_EEG_DF = os.sep.join([DATA_CACHE, f'{TODAY}-prediction-eeg-results.pkl'])
  840. if os.path.exists(results_eeg_df_name):
  841. print("Manuscript prediction results have been loaded.")
  842. elif os.path.exists(PREDICTION_RESULTS_EEG_DF):
  843. print(f"{PREDICTION_RESULTS_EEG_DF} prediction results have been loaded.")
  844. results_eeg_df = pd.read_pickle(PREDICTION_RESULTS_EEG_DF) # Prediction results
  845. else:
  846. for r in np.arange(n_runs):
  847. _results_df.append(get_results_df(data_v2, models_eeg[r], anes_data_keys, test_heldout[r], r))
  848. results_eeg_df = pd.concat(_results_df)
  849. pd.to_pickle(pd.DataFrame(results_eeg_df), PREDICTION_RESULTS_EEG_DF)
  853. # %% [markdown]
  854. # ### Facial features model
  856. # %%
  857. # def run_regression(ts_df, anes_data_keys, features, train_idx, alpha_range, models = None):
  858. def run_regression(ts_df, anes_data_keys, features, alpha_range):
  860. # specify feature parameters
  861. f_window = 1000 #6000 #In samples ; For moving_avg_filter the behavior features
  862. degree = 2 #3 #of the polynomial
  863. lag = 1
  865. # define empty lists to generate data for model
  866. measures = []
  867. time = []
  869. # data query for input features to the model
  870. meas_name_query = " | ".join(f"measurement_name == '{feature}'" for feature in features)
  872. indeces = []
  873. counter = 0
  874. for s_idx, dk in enumerate(anes_data_keys):
  875. d = dk[0]
  876. mouse = dk[1]
  878. anes_end = int(45 * 60 - ANES_START_OFFSETS[d, mouse])*100 # Trim at 45 minutes since anesthesia injection
  879. _measures = ts_df.query("mouse == @mouse & date == @d & condition == 'anes' & measurement_group != 'ephys' & "
  880. "(" + meas_name_query + ")").sort_values("measurement_name")["timeseries"].values
  882. # Filter noise
  883. for _ncol, _m in enumerate(_measures):
  884. _measures[_ncol] = moving_avg_filter(_m[:anes_end], f_window, stride = ds_factor[mouse, d])
  886. # Z-score
  887. _measures = z_score(np.stack(_measures, axis=1))
  888. _measures = get_features(_measures, degree, lag, interaction = True, initial_bias = True)
  889. _measures[:,0] = 1 #bias should be unaffected by zscoring and thus always 1
  891. # Save indeces in a list of arrays corresponding to each session
  892. indeces.append(np.arange(_measures.shape[0]) + counter)
  893. counter += _measures.shape[0]
  894. measures.append(_measures)
  896. # Define actual time
  897. ys = ts_df.query(
  898. "mouse == @mouse & "
  899. "condition == 'anes' & "
  900. "date == @d & "
  901. "measurement_name == @features[0]"
  902. )["timeseries"].values[0]
  903. _time = (np.arange(len(ys))/100 + ANES_START_OFFSETS[d, mouse])/ 60 #Note camera sample rate is 100
  904. _time = _time[:anes_end]
  905. _time = _time[f_window:-f_window]
  906. time.append(_time[::ds_factor[mouse, d]])
  908. measures = np.concatenate(measures, axis=0)
  909. time = np.concatenate(time, axis=0)
  910. cv = []
  911. # Annotate indeces to separate into validation and train sets
  912. for k in range(len(indeces)):
  913. train_indeces = np.concatenate([indeces[i] for i in range(len(indeces)) if i != k])
  914. test_indeces = indeces[k] #Leave one out cross-validation (a.k.a. test within the train set)
  915. cv.append((train_indeces, test_indeces))
  917. print(f"Running regression for {features}...")
  919. # fit Lasso model and explore alpha space
  920. mdl = Lasso(fit_intercept=True)
  921. param_grid = {'alpha': alpha_range} #key must be the same name that is used in Lasso documentation
  922. model = GridSearchCV(mdl, param_grid, cv = cv, verbose = 0, n_jobs = 6, scoring = "neg_root_mean_squared_error") # Explore hyper-parameter space (only for alpha in this case)
  923. model.fit(measures[:,1:], time)
  925. return model
  927. # Select facial features as input to the model
  928. # features = list(set(data_v2.query("measurement_group != 'ephys'")["measurement_name"]))
  929. features = ['eye-height-left', 'ear-angle-left', 'nose-bulge-volume']
  931. # Explore alpha space below:
  932. alpha_range = np.logspace(-1, 2, 100)
  933. n_runs = 220
  934. n_test = 3 #number of test sessions
  935. test_heldout =[]
  936. all_combos = list(itertools.combinations(range(len(anes_data_keys)), n_test))
  938. if os.path.exists(FF_MODEL_NAME):
  939. print("Manuscript model has been loaded.")
  940. elif os.path.exists(LASSO_MODEL_FF):
  941. print(f"{LASSO_MODEL_FF} has been loaded.")
  942. with open(LASSO_MODEL_FF, "rb") as fio:
  943. models_ff = pickle.load(fio) # Model results
  944. else:
  945. models_ff = []
  946. for r in np.arange(n_runs):
  947. test_heldout_idx = all_combos[r] #np.random.permutation(len(anes_data_keys))[:n_test]
  948. _test_heldout = [d for idx, d in enumerate(anes_data_keys) if idx in test_heldout_idx]
  949. print(test_heldout_idx, _test_heldout)
  950. test_heldout.append(_test_heldout)
  951. anes_data_keys_train = [d for d in anes_data_keys if d not in _test_heldout]
  952. models_ff.append(run_regression(data_v2, anes_data_keys_train, features, alpha_range))
  953. with open(LASSO_MODEL_FF, "wb") as fio:
  954. pickle.dump(models_ff, fio)
  958. # %%
  959. def get_results_df(data, model, anes_data_keys, test_heldout, features, ds_factor, n_run = 0):
  960. results = []
  961. # ANES_END = 10 * 60 * 100
  962. # These should match measurements used as input features
  963. meas_name_query = " | ".join(f"measurement_name == '{feature}'" for feature in features)
  964. for d in anes_data_keys:
  965. date = d[0]
  966. mouse = d[1]
  967. anes_end = int(45 * 60 - ANES_START_OFFSETS[date, mouse]) * 100# Trim at 45 minutes since anesthesia injection
  968. test_ts = data.query("date == @date & mouse == @mouse & condition == 'anes' & "
  969. "measurement_group != 'ephys' & "
  970. "(" + meas_name_query + ")").sort_values("measurement_name")["timeseries"].values
  972. # Filter noise
  973. f_window = 1000 #6000 # should be the same as the one used for the regression input features
  974. for _ncol, _m in enumerate(test_ts):
  975. test_ts[_ncol] = moving_avg_filter(_m[:anes_end], f_window, stride = ds_factor[mouse, date])
  977. # Z-score
  978. test_ts = z_score(np.stack(test_ts, axis=1))
  979. test_ts = get_features(test_ts, degree = 2, lags = 1, interaction = True, initial_bias = True)
  980. # test_ts = z_score(test_ts)
  981. test_ts[:,0] = 1
  983. # Define actual time
  984. ys = data.query(
  985. "mouse == @mouse & "
  986. "date == @date &"
  987. "condition == 'anes' & "
  988. "measurement_name == @features[0]"
  989. )["timeseries"].values[0]
  990. time = (np.arange(len(ys))/100 + ANES_START_OFFSETS[date, mouse])/ 60 #Note camera sample rate is 100
  991. time = time[:anes_end]
  992. time = time[f_window:-f_window]
  993. time = time[::ds_factor[mouse, date]]
  994. # predict time
  995. # m_idx = [0 if mouse == 'B47' else 1 if mouse == 'B48' else 2]
  996. # model = models[m_idx[0]]
  997. times_hat = model.predict(test_ts[:,1:])
  999. # Compute RMSE
  1000. rmse = np.sqrt(np.mean((times_hat - time) ** 2))
  1002. # save results
  1003. set_group = "Test" if (date, mouse, 'anes') in test_heldout else "Train"
  1004. results.append([n_run, mouse, date, set_group, "whole-face", times_hat, time, rmse])
  1006. results_df = pd.DataFrame(results, columns = ["run",
  1007. "mouse",
  1008. "date",
  1009. "set",
  1010. "features",
  1011. "predicted time",
  1012. "actual time",
  1013. "RMSE"])
  1014. return results_df
  1017. _results_df = []
  1018. PREDICTION_RESULTS_FF_DF = f'{TODAY}-prediction-ff-results.pkl'
  1021. if os.path.exists(results_df_name):
  1022. print("Manuscript model has been loaded.")
  1023. elif os.path.exists(PREDICTION_RESULTS_FF_DF):
  1024. print(f"{PREDICTION_RESULTS_FF_DF} has been loaded.")
  1025. results_df = pd.read_pickle(PREDICTION_RESULTS_FF_DF) # Prediction results
  1026. else:
  1027. for r in np.arange(n_runs):
  1028. _results_df.append(get_results_df(data_v2, models_ff[r], anes_data_keys, test_heldout[r], features, ds_factor, r))
  1029. results_df = pd.concat(_results_df)
  1030. pd.to_pickle(pd.DataFrame(results_df), PREDICTION_RESULTS_FF_DF)
  1033. # %% [markdown]
  1034. # ### figure
  1036. # %%
  1037. ex_mouse = 'B47'
  1038. ex_date = '20240912'
  1039. ex_run = [5, 7] #[face, eeg] #215
  1041. fig, ax = plt.subplots(nrows= 1,
  1042. ncols= 1,
  1043. figsize=(6,6),
  1044. dpi=200,
  1045. constrained_layout=True)
  1047. plot_predictions(results_eeg_df.query("run == @ex_run[1] & mouse == @ex_mouse & date == @ex_date")["actual time"],
  1048. results_eeg_df.query("run == @ex_run[1] & mouse == @ex_mouse & date == @ex_date")["predicted time"],
  1049. FILTER_WINDOW = 10,
  1050. color = EEG_CONTROL_CMAP[0],
  1051. label = 'eeg')
  1052. plot_predictions(results_df.query("run == @ex_run[0] & mouse == @ex_mouse & date == @ex_date")["actual time"],
  1053. results_df.query("run == @ex_run[0] & mouse == @ex_mouse & date == @ex_date")["predicted time"],
  1054. FILTER_WINDOW = 10,
  1055. color = CONTROL_CMAP[2],
  1056. label = 'facial features',
  1057. dash_line = False)
  1059. ax.set(title = f"{ex_mouse} ({ex_date})",
  1060. xlabel = "Actual time (min)", ylabel = "Predicted Time (min)",
  1061. xticks=[0, 10, 20, 30, 40, 50], yticks=[0, 10, 20, 30, 40, 50],
  1062. xlim = [0, 50], ylim = [0, 50])
  1063. ax.set_aspect("equal", "box")
  1064. sns.despine()
  1065. fname = f"3e-eeg-vs-ff-regression-ex-{ex_mouse}-{ex_date}"
  1066. fig.savefig(os.sep.join([OUTPUT_DIR, f"{fname}.svg"]), bbox_inches="tight")
  1068. # %% [markdown]
  1069. # ### statistics
  1071. # %%
  1072. # Save the mean RMSE results into .csv files
  1073. # mean_results_ff_df = test_mean_per_run(results_eeg_df, 220, "(subdelta, delta, theta)")
  1074. # mean_results_eeg_df = test_mean_per_run(results_df, 220, "whole-face)")
  1076. # mean_results_ff_df.to_csv(f'measurements-data-cache/CSVs/mean_results_ff_df_{TODAY}.csv', index=False)
  1077. # mean_results_eeg_df.to_csv(f'measurements-data-cache/CSVs/mean_results_eeg_df_{TODAY}.csv', index=False)
  1079. # %%
  1080. def cv_ttest_corrected(x, y, kfolds, nrepeats, ntrain, ntest):
  1081. diff = x - y
  1082. print(f"Length of the diff ({len(diff)}) should be equal to kfolds * nrepeats ({kfolds * nrepeats})")
  1083. v = np.sum((diff - np.mean(diff))**2) / (len(diff) - 1) # ~= np.std(diff, ddof=1)
  1084. tstat = np.mean(diff) / np.sqrt(v * (1/(kfolds * nrepeats) + ntest/ntrain)) # tstat with correction
  1085. pval = stats.t.sf(np.abs(tstat), nrepeats*kfolds - 1) # pvalue with correction
  1087. return tstat, pval
  1089. # %%
  1090. eeg_mean_results_df = test_mean_per_run(results_eeg_df, n_runs, "(subdelta, delta, theta)")
  1091. face_mean_results_df = test_mean_per_run(results_df, n_runs, "(subdelta, delta, theta)")
  1093. stat = cv_ttest_corrected(face_mean_results_df['mean_RMSE'].values,
  1094. eeg_mean_results_df['mean_RMSE'].values,
  1095. n_runs, 1, 9, 3)
  1096. print(f" tstat = {stat[0]}\n p-value = {stat[1]}")
  1098. # %% [markdown]
  1099. # ## Figure 3f
  1101. # %% [markdown]
  1102. # Root-mean-square error (RMSE) of time prediction where each dot represents the mean test error for one particular model trained on either Cheese3D (orange) or EEG (blue) features.
  1104. # %%
  1105. # Plot the results
  1106. fig, ax = plt.subplots(figsize=(3, 5))
  1108. # Facial features
  1109. n_runs = 220
  1110. sns.stripplot(test_mean_per_run(results_df, n_runs, "whole-face"),
  1111. x="features", y="mean_RMSE",
  1112. hue="features", #date #mouse #run
  1113. palette=CONTROL_CMAP_V2,
  1114. legend = False,
  1115. alpha=0.3, ax=ax)
  1116. sns.violinplot(test_mean_per_run(results_df, n_runs, "whole-face"),
  1117. x="features", y="mean_RMSE",
  1118. hue="features", #date #mouse #run
  1119. palette=CONTROL_CMAP_V2,
  1120. legend = False,
  1121. inner = None,
  1122. alpha=0.6, ax=ax)
  1125. # EEG power bands
  1126. sns.stripplot(test_mean_per_run(results_eeg_df, n_runs, "(subdelta, delta, theta)"),
  1127. x="features", y="mean_RMSE",
  1128. hue="features",
  1129. palette=EEG_CONTROL_CMAP,
  1130. legend = False,
  1131. alpha=0.3, ax=ax)
  1132. sns.violinplot(test_mean_per_run(results_eeg_df, n_runs, "(subdelta, delta, theta)"),
  1133. x="features", y="mean_RMSE",
  1134. hue="features",
  1135. palette=EEG_CONTROL_CMAP,
  1136. legend = False,
  1137. inner = None,
  1138. alpha=0.6, ax=ax)
  1141. ax.set_xlabel("Feature Set")
  1142. ax.set_ylabel("Root Mean Squared Error (min)")
  1143. ax.set_title("Prediction Error Across Feature Sets")
  1144. ax.set_xticks(list(range(len(ax.get_xticklabels()))),
  1145. [format_feature_ticks(f) for f in ax.get_xticklabels()])
  1146. ax.tick_params(axis="x", bottom=False)
  1147. ax.set_ylim(0, None)
  1148. ax.set_yticks(np.arange(0,16,5))
  1149. ax.legend(frameon=False, bbox_to_anchor=(1.0, 0.3))
  1150. sns.despine(fig)
  1152. fname = f"3f-eeg-vs-face-regression-results-summary-violin"
  1153. fig.savefig(os.sep.join([OUTPUT_DIR, f"{fname}.svg"]), bbox_inches="tight")
  1155. # %% [markdown]
  1156. # ## Figure 3g
  1158. # %% [markdown]
  1159. # (notebook: 2024-12-05-long-anes-three-way-predictions-eeg)
  1160. #
  1161. # Overlaid z-scored theta, delta, and sub-delta frequency band power traces from 12 sessions (same as in (c);
  1162. # vertical scale bars indicate one standard deviation)
  1164. # %%
  1165. eeg_window = 8
  1166. spacer = 8 #5
  1167. features = ["subdelta (<1Hz)", "delta [1-4Hz]", "theta [5-10Hz]"]
  1168. scale = 1
  1169. fig, axs = plt.subplots(nrows=len(features),
  1170. ncols=1,
  1171. figsize=(5, 8),
  1172. dpi=150,
  1173. constrained_layout=True)
  1175. meas_summary = {feat: [] for feat in features}
  1176. for sess in anes_data_keys:
  1178. mouse = sess[1]
  1179. date = sess[0]
  1180. eegpower = data_v2.query("mouse == @mouse & date == @date & condition == 'anes' & measurement_name == 'eeg-power'")["timeseries"]
  1181. fs, ts, power = eegpower.values[0]
  1182. bands = compute_power_bands(fs, power) #type(_bands) = tuple
  1183. anes_start = np.max(list(ANES_START_OFFSETS_EEG.values())) - ANES_START_OFFSETS_EEG[date, mouse] #in seconds - Take the same start for all of them
  1184. anes_start = int(np.floor(anes_start) * eeg_power_sr[mouse, date])
  1185. anes_end = 540 #int(np.floor((45 * 60 - anes_start) * eeg_power_sr[mouse, date])) # Trim at 45 minutes since anesthesia injection
  1187. for i_b, _b in enumerate(bands):
  1188. # FILTER AND Z SCORED DATA PER MEASUREMENT
  1189. if i_b > 2:
  1190. continue
  1191. shortened_b = _b[anes_start:anes_start + anes_end]
  1192. time = np.arange(0, len(shortened_b))/eeg_power_sr[mouse, date] + np.max(list(ANES_START_OFFSETS_EEG.values()))
  1193. time = (time[eeg_window:-eeg_window])/60
  1194. # print(time[-1], len(time))
  1195. ys = z_score(moving_avg_filter(shortened_b, window = eeg_window, stride = 1))
  1196. sns.lineplot(x = time, y = ys,
  1197. color = EEG_CONTROL_CMAP[0], linewidth = 0.5, alpha = 0.5,
  1198. ax = axs[i_b])
  1199. meas_summary[features[i_b]].append(ys)
  1200. axs[i_b].sharex(axs[0])
  1201. axs[i_b].set_ylabel(('\n').join(features[i_b].split(' ')), rotation = 0)
  1202. axs[i_b].yaxis.set_label_coords(-0.2,0.5)
  1203. if i_b != len(features)-1:
  1204. axs[i_b].tick_params("x", bottom=False, labelbottom=False)
  1205. sns.despine(ax=axs[i_b], bottom=True)
  1206. # axs[col].plot(time[f_window:-f_window], ys - spacer*spacing, color = MEASUREMENT_CMAP[l])
  1207. # plt.plot(time, ys - spacer*spacing, color = MEASUREMENT_CMAP[l])
  1208. for row, feat in enumerate(features):
  1209. sns.lineplot(x = time, y = np.mean(meas_summary[feat], axis=0),
  1210. color = EEG_CONTROL_CMAP[0], linewidth = 2.5, ax = axs[row])
  1212. sns.despine(ax = axs[2])
  1214. axs[len(features)-1].set_xlabel('Time (s)')
  1215. fig.tight_layout()
  1216. # axs[0].set_yticks(np.arange(-spacer,-(spacer*spacing+1), -spacer)) #Negative signs to invert the order of measurements top-bottom
  1217. # axs[0].set_yticklabels(y_labels);
  1219. file_name = "eeg-band-feat-with-mean"
  1220. # save_figure("rev_figs/", f"2024-{file_name}", fig, formats=["svg"])
  1221. # fig.savefig(f'{file_name}_ZSCORED.pdf')
  1223. # %% [markdown]
  1224. # ## Figure 3a
  1226. # %% [markdown]
  1227. # (notebook: 2023-long-anes-measurements.ipynb)
  1229. # %% [markdown]
  1230. # Example facial movement raster plot during anesthesia with concurrent EEG recording (each vertical line cor- responds to movement above the 99.9-th percentile jitter threshold as shown in Supplementary Figure 4c for a given 10 ms time window).
  1232. # %% [markdown]
  1233. # ### load data and functions from 2023 cohort
  1235. # %%
  1236. # For faster development, limit data to X mice
  1237. # Set to a high number to render plots for every mouse
  1238. MAX_NUM_MICE = 100
  1240. # Plot every row
  1241. PLOT_SAMPLE = 100
  1243. # Smooth data for visualization
  1244. SMOOTH_WINDOW_SIZE = 25
  1246. # We didn't measure the amount of time between the end of the awake video and the anesthesia injection
  1247. # This adds a constant offset of 5 minutes for every mouse
  1248. AWAKE_END_OFFSET = 5*60
  1250. # Start of anesthesia video after injection (in seconds)
  1251. ANES_START_OFFSETS = {
  1252. 'B6': 120,
  1253. 'B8': 180,
  1254. 'B15': 60,
  1255. 'B20': 31,
  1256. 'B26': 53,
  1257. 'B33': 30,
  1258. }
  1260. STILL_PERIODS = {
  1261. "B6": (38.50, 44.05),
  1262. "B8": (34.30, 39.55),
  1263. "B15": (27.20, 34.15),
  1264. "B20": (34.55, 41.00),
  1265. "B26": (58.20, 63.28),
  1266. "B33": (44.30, 50.05) #This mouse was commented out in the original notebook (2024-04-05-long-anes-jitter.ipynb)
  1267. # "C3": (0, None)
  1268. }
  1270. _, MEASUREMENT_CMAP = measurements_cmap()
  1272. DATA_CACHE = "measurements-data-cache-2023"
  1273. os.makedirs(DATA_CACHE, exist_ok=True)
  1274. COORD_DATA_CACHE = os.sep.join([DATA_CACHE, "long-anes-coords.pkl"])
  1275. MEASURE_DATA_CACHE = os.sep.join([DATA_CACHE, "long-anes-meas.pkl"])
  1276. JITTER_DATA_CACHE = os.sep.join([DATA_CACHE, "long-anes-jitter-meas.pkl"])
  1278. ANIPOSE_BASE = 'anipose-projects/20231013-long-anes-rig2'
  1279. COORDINATE_PATHS = {}
  1280. key_cols = ('mouse', 'source', 'condition')
  1281. for p in Path(ANIPOSE_BASE).glob('*/pose-3d/*.csv'):
  1282. mouse = p.name.split('_')[1]
  1283. source = 'rig2'
  1284. condition = 'awake' if 'awake' in p.name else 'anes'
  1285. COORDINATE_PATHS[(mouse, source, condition)] = p
  1286. data_keys = list(COORDINATE_PATHS.keys())
  1287. data_keys
  1289. # %%
  1290. def build_measjitter_df(meas_df, periods):
  1291. # select only the subset of rows that match the mouse/source/condition pairs
  1292. queries = [(meas_df["mouse"] == mouse) &
  1293. (meas_df["source"] == "rig2") &
  1294. ((meas_df["condition"] == "anes") | (meas_df["condition"] == "dead"))
  1295. for mouse in periods.keys()]
  1296. sub_df = meas_df[reduce(lambda x, y: x | y, queries)]
  1297. jitter_df = sub_df[[*key_cols,
  1298. "measurement_group",
  1299. "measurement_name",
  1300. "measurement_type",
  1301. "timeseries"]]
  1302. jitter_df = jitter_df.rename(columns={"measurement_group": "region"})
  1303. jitter_df.loc[(jitter_df["region"] == 'eye') &
  1304. (jitter_df["measurement_name"].str.contains('left')), "region"] = "eye(left)"
  1305. jitter_df.loc[(jitter_df["region"] == 'eye') &
  1306. (jitter_df["measurement_name"].str.contains('right')), "region"] = "eye(right)"
  1307. jitter_df.loc[(jitter_df["region"] == 'ear') &
  1308. (jitter_df["measurement_name"].str.contains('left')), "region"] = "ear(left)"
  1309. jitter_df.loc[(jitter_df["region"] == 'ear') &
  1310. (jitter_df["measurement_name"].str.contains('right')), "region"] = "ear(right)"
  1311. jitter_df.loc[(jitter_df["region"] == "cheek"), "region"] = "whisker pad"
  1312. jitter_df.loc[(jitter_df["measurement_name"] == "cheek-bulge-volume"), "measurement_name"] = "cheek-bulge-volume"
  1313. jitter_df.loc[(jitter_df["measurement_name"] == "nose-bulge-volume"), "measurement_name"] = "nose-bulge-volume"
  1314. for mouse, period in periods.items():
  1315. start = round(period[0] * 60 * 100) if period[0] is not None else None
  1316. end = round(period[1] * 60 * 100) if period[1] is not None else None
  1317. idx = jitter_df["mouse"] == mouse
  1318. jitter_df.loc[idx, "timeseries"] = jitter_df.loc[idx, "timeseries"].apply(
  1319. lambda x: x[start:end]
  1320. )
  1321. jitter_df["deviations"] = jitter_df.groupby("measurement_name")["timeseries"].transform(
  1322. lambda x: x.apply(lambda xi: xi - xi.mean())
  1323. )
  1324. jitter_df["timeseries_stddev"] = jitter_df.groupby("measurement_name")["deviations"].transform(
  1325. lambda x: x.apply(np.std)
  1326. )
  1327. jitter_df["min_deviation"] = jitter_df.groupby("measurement_name")["deviations"].transform(
  1328. lambda x: x.apply(np.min)
  1329. )
  1330. jitter_df["max_deviation"] = jitter_df.groupby("measurement_name")["deviations"].transform(
  1331. lambda x: x.apply(np.max)
  1332. )
  1333. jitter_df["velocities"] = jitter_df.groupby("measurement_name")["timeseries"].transform(
  1334. lambda x: x.apply(lambda xi: np.abs(np.diff(xi))) * 100
  1335. )
  1336. jitter_df["velocity_mean"] = jitter_df.groupby("measurement_name")["velocities"].transform(
  1337. lambda x: x.apply(np.mean)
  1338. )
  1339. jitter_df["velocity_thresh"] = jitter_df.groupby("measurement_name")["velocities"].transform(
  1340. lambda x: x.apply(lambda xi: np.percentile(xi, 99.9))
  1341. )
  1342. jitter_df["min_velocity"] = jitter_df.groupby("measurement_name")["velocities"].transform(
  1343. lambda x: x.apply(np.min)
  1344. )
  1345. jitter_df["max_velocity"] = jitter_df.groupby("measurement_name")["velocities"].transform(
  1346. lambda x: x.apply(np.max)
  1347. )
  1348. jitter_df["deviations_au"] = jitter_df.groupby("measurement_name")["timeseries"].transform(
  1349. lambda x: x.apply(lambda xi: (xi - xi.mean()) / xi.mean())
  1350. )
  1351. jitter_df["timeseries_stddev_au"] = jitter_df.groupby("measurement_name")["deviations_au"].transform(
  1352. lambda x: x.apply(np.std)
  1353. )
  1354. jitter_df["min_deviation_au"] = jitter_df.groupby("measurement_name")["deviations_au"].transform(
  1355. lambda x: x.apply(np.min)
  1356. )
  1357. jitter_df["max_deviation_au"] = jitter_df.groupby("measurement_name")["deviations_au"].transform(
  1358. lambda x: x.apply(np.max)
  1359. )
  1360. jitter_df["velocities_au"] = jitter_df.groupby("measurement_name")["timeseries"].transform(
  1361. lambda x: x.apply(lambda xi: np.abs(np.diff(xi)) / np.mean(np.abs(np.diff(xi))))
  1362. )
  1363. jitter_df["velocity_mean_au"] = jitter_df.groupby("measurement_name")["velocities_au"].transform(
  1364. lambda x: x.apply(np.mean)
  1365. )
  1366. jitter_df["min_velocity_au"] = jitter_df.groupby("measurement_name")["velocities_au"].transform(
  1367. lambda x: x.apply(np.min)
  1368. )
  1369. jitter_df["max_velocity_au"] = jitter_df.groupby("measurement_name")["velocities_au"].transform(
  1370. lambda x: x.apply(np.max)
  1371. )
  1373. return jitter_df
  1375. # %%
  1376. from scipy.signal import medfilt
  1378. if not os.path.exists(COORD_DATA_CACHE):
  1379. print(f"Pre-filtering coordinate data and storing in {COORD_DATA_CACHE}...")
  1380. coord_data = {k: read_3d_data(v.parent.parent.as_posix(),
  1381. filter_func=medfilt,
  1382. filter_kwargs={'kernel_size': (SMOOTH_WINDOW_SIZE, 1)})
  1383. for k, v in COORDINATE_PATHS.items()}
  1384. with open(COORD_DATA_CACHE, 'wb') as dict_pkl:
  1385. pickle.dump(coord_data, dict_pkl)
  1386. else:
  1387. with open(COORD_DATA_CACHE, 'rb') as dict_pkl:
  1388. coord_data = pickle.load(dict_pkl)
  1390. if not os.path.exists(MEASURE_DATA_CACHE):
  1391. print(f"Pre-computing measurements data and storing in {MEASURE_DATA_CACHE}...")
  1392. meas_df = compute_measurements_df(coord_data)
  1393. meas_df.to_pickle(MEASURE_DATA_CACHE)
  1394. else:
  1395. meas_df = pd.read_pickle(MEASURE_DATA_CACHE)
  1397. if os.path.exists(JITTER_DATA_CACHE):
  1398. jitter_df = pd.read_pickle((JITTER_DATA_CACHE))
  1399. else:
  1400. jitter_df = build_measjitter_df(meas_df, STILL_PERIODS)
  1401. # raise RuntimeError("Jitter results not pre-computed. Run jitter analysis notebook.")
  1403. # %%
  1404. awake_end = {
  1405. mouse: len(df)
  1406. for (mouse, _, condition), df in coord_data.items()
  1407. if condition == 'awake'
  1408. }
  1409. awake_end
  1411. # %%
  1412. from collections import defaultdict
  1414. mice = list(meas_df["mouse"].unique())
  1415. measurement_names = list(meas_df['measurement_name'].unique())
  1416. measurement_group_names = list(meas_df['measurement_group'].unique())
  1417. measurement_groups = defaultdict(set)
  1418. for measurement_name in measurement_names:
  1419. sub_df = meas_df.query("measurement_name == @measurement_name")
  1420. measurement_groups[sub_df['measurement_group'].iloc[0]].add(sub_df['measurement_name'].iloc[0])
  1421. measurement_groups
  1423. # %%
  1424. from labutils.utils import unzip
  1426. def build_timeseries_df(meas_df, data_keys):
  1427. # select only the subset of rows that match the mouse/source/condition pairs
  1428. queries = [(meas_df["mouse"] == mouse) &
  1429. (meas_df["source"] == source) &
  1430. (meas_df["condition"] == condition)
  1431. for mouse, source, condition in data_keys]
  1432. sub_df = meas_df[reduce(lambda x, y: x | y, queries)]
  1433. # compute velocities and pad
  1434. def _process_timeseries(row):
  1435. v = np.abs(np.diff(row["timeseries"]))
  1436. if row["condition"] == "awake":
  1437. p = np.concatenate([row["timeseries"],
  1438. np.zeros(100 * AWAKE_END_OFFSET)])
  1439. v = np.concatenate([[0], v, np.zeros(100 * AWAKE_END_OFFSET)])
  1440. t = np.arange(-len(p), 0)
  1441. else:
  1442. p = np.concatenate([np.zeros(100 * ANES_START_OFFSETS[row["mouse"]]),
  1443. row["timeseries"]])
  1444. v = np.concatenate([np.zeros(100 * ANES_START_OFFSETS[row["mouse"]] + 1), v])
  1445. t = np.arange(len(p))
  1447. return p, v, t
  1448. positions, velocities, frames = unzip(sub_df.apply(_process_timeseries, axis=1).values)
  1449. # create new dataframe for velocities spread out over time
  1450. ts_df = sub_df[[*key_cols, "measurement_name"]].copy()
  1451. ts_df = ts_df.assign(position=positions, velocity=velocities, frames=frames)
  1452. ts_df = ts_df.explode(["position", "velocity", "frames"])
  1453. ts_df = ts_df.pivot(index=[*key_cols, "frames"],
  1454. columns="measurement_name",
  1455. values=["position", "velocity"])
  1456. ts_df.columns = ["/".join(multicol).strip() for multicol in ts_df.columns.values]
  1457. ts_df.reset_index(inplace=True)
  1458. ts_df = ts_df.assign(**{"time (s)": ts_df["frames"] / 100,
  1459. "time (min)": ts_df["frames"] / 100 / 60})
  1460. for name in measurement_names:
  1461. ts_df[f"velocity/{name}-std"] = ts_df.groupby("mouse")[f"velocity/{name}"].transform(
  1462. lambda x: x / x.std()
  1463. )
  1464. ts_df.sort_values(["mouse", "frames"], inplace=True)
  1466. return ts_df
  1468. ts_df = build_timeseries_df(meas_df, data_keys)
  1469. ts_df
  1471. # %% [markdown]
  1472. # ### figure
  1474. # %%
  1475. PLOT_SAMPLE = 1
  1477. time_col = "time (min)"
  1478. plot_meas = ["velocity/ear-angle-left",
  1479. "velocity/ear-angle-right",
  1480. "velocity/eye-area-left",
  1481. "velocity/eye-area-right",
  1482. "velocity/mouth-area",
  1483. "velocity/cheek-bulge-volume",
  1484. "velocity/nose-bulge-volume"]
  1485. cols = [time_col] + plot_meas
  1486. column_order = [time_col] + plot_meas
  1487. fig, axs = plt.subplots(nrows=len(mice),
  1488. # sharex=True,
  1489. figsize=(10, 3 * len(mice)))
  1491. for mouse, ax in zip(filter(lambda x: x != "B33", mice), axs):
  1492. # get plotting subset
  1493. sub_df = ts_df.query('mouse == @mouse')[cols].copy()
  1494. sub_df = sub_df[column_order]
  1495. sub_df.set_index(time_col, inplace=True)
  1496. sub_df = sub_df.iloc[::PLOT_SAMPLE, :]
  1498. # ax = axs[ax_num]
  1499. for i, column in enumerate(plot_meas):
  1500. meas_name = column.split("/")[-1]
  1501. baseline = jitter_df.query("mouse == @mouse & "
  1502. "measurement_name == @meas_name")["velocities"]
  1503. baseline = baseline.values[0] / 100 # convert mm / s -> mm / frame
  1504. thresh = np.percentile(baseline, 99.9)
  1505. # thresh = np.std(baseline)
  1506. print(mouse, meas_name, "threshold =", thresh)
  1507. ticks = sub_df.index[sub_df[column] > thresh]
  1508. # Plot vertical lines for threshold crossing events
  1509. if len(ticks) > 0:
  1510. ax.vlines(ticks, i, i - 1,
  1511. colors=MEASUREMENT_CMAP[meas_name],
  1512. alpha=0.01,
  1513. linewidths=0.5)
  1515. ax.set_title(mouse)
  1516. ax.tick_params(axis="y", left=False)
  1517. ax.set_yticklabels([])
  1518. ax.get_xaxis().set_ticks_position('bottom')
  1519. sns.despine(fig, left=True, top=True, right=True)
  1520. fig.tight_layout()
  1522. fname = f"3a-long-anes-tick-rasters-allmice"
  1523. # fig.savefig(os.sep.join([OUTPUT_DIR, f"{fname}.png"]), bbox_inches="tight")
  1524. # save_figure(OUTPUT_DIR, "long-anes-panel-e-exemplar_ticks", fig, formats=["png"])
  1526. # %% [markdown]
  1527. # ## Figure 3b
  1529. # %% [markdown]
  1530. # (notebook: 2023-long-anes-measurements.ipynb)
  1532. # %% [markdown]
  1533. # Zoom in of movement raster plot from (a) to show the early moments of movement recovery following anesthesia.
  1535. # %%
  1536. fig, axs = plt.subplots(1, len(mice), figsize=(2 * len(mice), 2))
  1538. wakeup_window = {
  1539. "B15": (48, 52),
  1540. "B20": (62, 66),
  1541. "B26": (101, 105),
  1542. "B6": (48, 52),
  1543. "B8": (48, 52)
  1544. }
  1546. for mouse, ax in zip(filter(lambda x: x != "B33", mice), axs):
  1547. sub_df = ts_df.query('mouse == @mouse')[cols].copy()
  1548. sub_df = sub_df[column_order]
  1549. sub_df.set_index(time_col, inplace=True)
  1550. sub_df = sub_df.iloc[::PLOT_SAMPLE, :]
  1552. for i, column in enumerate(plot_meas):
  1553. meas_name = column.split("/")[-1]
  1554. baseline = jitter_df.query("mouse == @mouse & "
  1555. "measurement_name == @meas_name")["velocities"]
  1556. baseline = baseline.values[0] / 100 # convert mm / s -> mm / frame
  1557. thresh = np.percentile(baseline, 99.9)
  1558. # thresh = np.std(baseline)
  1559. print(mouse, meas_name, "threshold =", thresh)
  1560. ticks = sub_df.index[sub_df[column] > thresh]
  1561. # Plot vertical lines for threshold crossing events
  1562. if len(ticks) > 0:
  1563. ax.vlines(ticks, i, i - 1,
  1564. colors=MEASUREMENT_CMAP[meas_name],
  1565. alpha=0.1,
  1566. linewidths=0.5)
  1568. ax.set_xlim(*wakeup_window[mouse])
  1569. ax.set_title(mouse)
  1570. ax.tick_params(axis="y", left=False)
  1571. ax.set_yticklabels([])
  1572. ax.get_xaxis().set_ticks_position('bottom')
  1573. sns.despine(fig, left=True, top=True, right=True)
  1574. fig.tight_layout()
  1576. fname = f"3b-long-anes-tick-rasters-allmice-zoom"
  1577. # fig.savefig(os.sep.join([OUTPUT_DIR, f"{fname}.png"]), bbox_inches="tight")
  1578. # save_figure(OUTPUT_DIR, "long-anes-panel-f-exemplar_ticks_zoom", fig, formats=["png"])

fig3-part1-cheese3d-general-anesthesia-eeg.ipynb at commit 07ef417, under MIT · at the source

Overview

Authors: Kyle Daruwalla1, Irene Nozal Martin1,2, Linghua Zhang1, Diana Naglič1, Andrew Frankel1, Catherine Rasgaitis1, Rubin Zhao1, Xinyan Zhang1, Zainab Ahmad1, Jeremy C Borniger1, Xun Helen Hou1
ORCID iDs: Irene Nozal Martin, Linghua Zhang, Catherine Rasgaitis, Jeremy C Borniger, Xun Helen Hou
  1. Cold Spring Harbor Laboratory, Cold Spring Harbor, NY USA
  2. Dept. of Neuroscience, Stony Brook University, Stony Brook, NY USA
Institutions: Cold Spring Harbor Laboratory (United States); Stony Brook University (United States)
Journal: Nature neuroscience, volume 29, issue 6, pages 1510-1521
Dates: received 3 February 2025; accepted 6 March 2026; published online 27 April 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41593-026-02262-8 · PMID 42045464 · PMCID PMC13246446 · OpenAlex W7155908500
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), systems (subfield)
Methods: Preprocessing, Statistics, Machine learning, Spectral & time-frequency, fMRI & imaging, Single-unit activity, calcium imaging, Smoothing, state filtering, decompositions, Physiology & signal measures, Connectivity
Keywords: Behavioural methods, Image processing, Hardware and infrastructure, Motor control, Neurophysiology
MeSH: Face*, Facial Expression*, Imaging, Three-Dimensional*, Animals, Male, Mice, Mice, Inbred C57BL, Motion Capture, Movement (* major topic)
Topic: Face Recognition and Perception (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Fulbright Association; Brain and Behavior Research Foundation (Young investigator Award #29769); NIBIB NIH HHS (T32 EB003383); Schmidt Foundation, CSHL NeuroAI Scholarship
Citations: cited by 4 papers (Europe PMC); 63 references in the paper

Abstract

Facial expressions and movements, from a subtle and ephemeral grimace to vigorous and rapid chewing, offer direct insights into the moment-to-moment changes of neural and physiological processes. Mice, with discernible facial responses and evolutionarily conserved mammalian facial movement control circuits, provide an ideal model in which to unravel the link between facial movement and underlying states. However, existing frameworks lack the spatial or temporal resolution to sensitively track all movements of the mouse face because of its small and conical form factor. We introduce Cheese3D, a computer vision system that captures high-speed 3D motion of the entire mouse face (including ears, eyes, whisker pad and jaw, covering both sides of the face), using a calibrated six-camera array. The interpretable framework extracts dynamics of anatomically meaningful 3D facial features in absolute world units at sub-mm precision. The precise face-wide motion data generated by Cheese3D provides clear insights, as shown by proof-of-principle experiments predicting anesthetic depth from changing facial patterns, inferring tooth and muscle anatomy from fast ingestion motions across the entire face, measuring minute differences in movements evoked by brainstem stimulation and relating neural activity to spontaneous facial movements, including expressive features only measurable in 3D (for example, angles of ear motion). Cheese3D can serve as a discovery tool that renders subtle mouse facial movements as a highly interpretable readout of otherwise hidden processes.

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 19 matches between paragraphs and lines of code.

Hou-Lab-CSHL/cheese3d

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 07ef417d37166ad1a272767340827ebe11bd1fe2, 7 July 2026
Languages: Python (27), Jupyter (14), Shell (1)
Size: 121 files, 42 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (pixi.lock, pixi.toml, public-env/pixi.lock, public-env/pixi.toml, packages/cheese3d-annotator/pyproject.toml, packages/cheese3d/pyproject.toml), tests, continuous integration, documentation, 14 notebooks
Not found: CITATION.cff
Tools: NumPy (20 files), pandas (16 files), Matplotlib (10 files), seaborn (7 files), SciPy (6 files), napari (5 files), OpenCV (5 files), scikit-learn (3 files), DeepLabCut (1 file), h5py (1 file), imageio (1 file), JAX (1 file), Open Ephys analysis tools (1 file), scikit-image (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
35 files
At the source: github.com/Hou-Lab-CSHL/cheese3d

Zenodo 18573618

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Languages: Jupyter (14)
Size: 15 files, 14 scripts
Software Heritage: not checked
Found in: “Code availability”
Holds: README, 14 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (10 files), NumPy (10 files), pandas (10 files), seaborn (10 files), SciPy (6 files), scikit-learn (4 files), h5py (2 files), JAX (1 file), Kilosort (1 file), OpenCV (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
11 files
At the source:

Code availability

Cheese3D software and example datasets are available at https://github.com/Hou-Lab-CSHL/cheese3d. Additionally, a copy of the analysis code is available on Zenodo at 10.5281/zenodo.18573618 (ref. 63).

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:

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

Datasets cited

Data availability

All raw data used for analyses in the paper are publicly available on Zenodo at 10.5281/zenodo.18508087 (ref. 62). 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, 11 authors, 5 keywords, 9 MeSH terms, 4 funders, 53 references.

Cite

This paper

Daruwalla, K., Nozal Martin, I., Zhang, L., Naglič, D., Frankel, A., Rasgaitis, C., Zhao, R., Zhang, X., Ahmad, Z., Borniger, J. C., & Hou, X. H. (2026). Cheese3D enables sensitive detection and analysis of whole-face movement in mice. Nature neuroscience, 29(6), 1510-1521. https://doi.org/10.1038/s41593-026-02262-8

BibTeX

@article{daruwalla2026cheese3d,
author = {Daruwalla, Kyle and Nozal Martin, Irene and Zhang, Linghua and Naglič, Diana and Frankel, Andrew and Rasgaitis, Catherine and Zhao, Rubin and Zhang, Xinyan and Ahmad, Zainab and Borniger, Jeremy C and Hou, Xun Helen},
title = {{Cheese3D enables sensitive detection and analysis of whole-face movement in mice}},
journal = {Nature neuroscience},
year = {2026},
month = apr,
volume = {29},
number = {6},
pages = {1510--1521},
publisher = {Nature Portfolio},
issn = {1097-6256},
doi = {10.1038/s41593-026-02262-8},
url = {https://doi.org/10.1038/s41593-026-02262-8},
pmid = {42045464},
pmcid = {PMC13246446}
}

RIS

TY - JOUR
AU - Daruwalla, Kyle
AU - Nozal Martin, Irene
AU - Zhang, Linghua
AU - Naglič, Diana
AU - Frankel, Andrew
AU - Rasgaitis, Catherine
AU - Zhao, Rubin
AU - Zhang, Xinyan
AU - Ahmad, Zainab
AU - Borniger, Jeremy C
AU - Hou, Xun Helen
TI - Cheese3D enables sensitive detection and analysis of whole-face movement in mice
T2 - Nature neuroscience
J2 - Nat Neurosci
PY - 2026
DA - 2026/04/27
VL - 29
IS - 6
SP - 1510
EP - 1521
SN - 1097-6256
PB - Nature Portfolio
DO - 10.1038/s41593-026-02262-8
UR - https://doi.org/10.1038/s41593-026-02262-8
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41593-026-02262-8",
"type": "article-journal",
"title": "Cheese3D enables sensitive detection and analysis of whole-face movement in mice",
"container-title": "Nature neuroscience",
"author": [
{
"family": "Daruwalla",
"given": "Kyle"
},
{
"family": "Nozal Martin",
"given": "Irene"
},
{
"family": "Zhang",
"given": "Linghua"
},
{
"family": "Naglič",
"given": "Diana"
},
{
"family": "Frankel",
"given": "Andrew"
},
{
"family": "Rasgaitis",
"given": "Catherine"
},
{
"family": "Zhao",
"given": "Rubin"
},
{
"family": "Zhang",
"given": "Xinyan"
},
{
"family": "Ahmad",
"given": "Zainab"
},
{
"family": "Borniger",
"given": "Jeremy C"
},
{
"family": "Hou",
"given": "Xun Helen"
}
],
"container-title-short": "Nat Neurosci",
"volume": "29",
"issue": "6",
"page": "1510-1521",
"DOI": "10.1038/s41593-026-02262-8",
"PMID": "42045464",
"PMCID": "PMC13246446",
"ISSN": "1097-6256",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41593-026-02262-8",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
27
]
]
}
}

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.celrep.2026.117420 [code]
Neural population dynamics of direct electrical stimulation of neocortex.
Journal: Cell reports
In common: Open Ephys analysis tools, napari, JAX, 9 other tools, systems, mouse
[2] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: DeepLabCut, Kilosort, OpenCV, 8 other tools, systems, mouse, 2 references
[3] doi:10.7554/elife.109717 [code]
Retrosplenial cortex enables context-dependent goal-directed sensorimotor transformation.
Journal: eLife
In common: Kilosort, imageio, OpenCV, 8 other tools, systems, mouse, 1 reference
[4] doi: [code]
Real-time closed-loop feedback system for mouse mesoscale cortical signal and movement control
Journal: eLife
In common: napari, imageio, OpenCV, 8 other tools, mouse, 1 reference
[5] doi:10.1038/s41467-026-72152-x [code]
Centralized brain networks controlling antennal grooming coordination.
Journal: Nature communications
In common: DeepLabCut, OpenCV, h5py, 6 other tools, systems, 2 references
[6] doi:10.1038/s41586-026-10679-1 [code]
Cortical development dynamics across autism spectrum disorder mouse models.
Journal: Nature
In common: imageio, OpenCV, scikit-image, 7 other tools, mouse, 1 reference
[7] doi:10.1038/s41467-026-72057-9 [code]
Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.
Journal: Nature communications
In common: OpenCV, scikit-image, h5py, 6 other tools, systems, 2 references
[8] doi:10.1016/j.crmeth.2026.101421 [code]
EthoPy provides an accessible platform for reproducible behavioral neuroscience.
Journal: Cell reports methods
In common: imageio, OpenCV, h5py, 6 other tools, mouse, 2 references
[9] doi:10.3389/fendo.2026.1828487 [code]
Castration-induced nigrostriatal deficits are linked to reduced TrkB and loss of mature spines in the dorsal striatum.
Journal: Frontiers in endocrinology
In common: imageio, OpenCV, scikit-image, 7 other tools, mouse, 1 reference
[10] 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: napari, imageio, OpenCV, 7 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.