OSCR

Fear conditioning biases olfactory sensory neuron frequencies across generations.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

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 · 700 lines · 24 KB · MIT

  1. # %%
  2. # import necessary packages
  3. import os
  4. import numpy as np
  5. import pandas as pd
  6. import seaborn as sns
  7. import matplotlib.pyplot as plt
  8. from scipy import stats
  9. from statsmodels.stats.multitest import multipletests
  10. from itertools import combinations
  11. from matplotlib import rc
  12. rc('font',**{'family':'sans-serif','sans-serif':['Arial']})
  14. import warnings
  15. warnings.filterwarnings('ignore')
  17. # %%
  18. # specify pathnames of directories and moseq_df/stats_df csvs generated by kpMoSeq
  20. project_dir= '' # the full path to the project directory
  21. model_name='' # name of model to analyze (e.g. something like `2023_05_23-15_19_03`)
  22. fig_dir = '' # name of folder you want figures to be saved in
  23. stats_dir = project_dir+model_name+'/stats/' # where you want stats files to be saved e.g. significant syllables, p-vals of pairwise comparisons
  24. moseq_df_filename = 'moseq_df.csv' # csv filename of saved moseq_df
  25. stats_df_filename = 'p_up_p_stats_df.csv' # csv filename of saved stats_df
  26. sig_syllables_filename = 'p_up_p_sig_syllables_freq.csv' # csv filename of saved sig_syllables df, if saved -- if not, can generate below
  28. moseq_df = pd.read_csv(project_dir+model_name+'/moseq_df/'+moseq_df_filename) # read in moseq_df
  29. stats_df = pd.read_csv(project_dir+model_name+'/stats_df/'+stats_df_filename) # read in stats_df
  31. # %%
  32. # define groups for comparison (group_1, group_2), statistic of interest (frequency or duration), and p-val
  33. group_1 = 'f1_up'
  34. group_2 = 'f1_p'
  35. stat = 'frequency'
  36. thresh = 0.05
  38. # %%
  39. # define functions needed for analysis (from kpMoSeq code)
  41. def run_kruskal(
  42. stats_df,
  43. statistic="frequency",
  44. n_perm=10000,
  45. seed=42,
  46. thresh=0.05,
  47. mc_method="fdr_bh",
  48. ):
  49. """Run Kruskal-Wallis test on syllable usage data.
  51. Parameters
  52. ----------
  53. stats_df : pandas.DataFrame
  54. DataFrame containing syllable usage data.
  55. statistic : str, optional
  56. Statistic to use for KW test, by default 'frequency'
  57. n_perm : int, optional
  58. Number of permutations to run, by default 10000
  59. seed : int, optional
  60. Random seed, by default 42
  61. thresh : float, optional
  62. Alpha threshold to consider syllable significant, by default 0.05
  63. mc_method : str, optional
  64. Multiple Corrections method to use, by default "fdr_bh"
  66. Returns
  67. -------
  68. df_k_real : pandas.DataFrame
  69. DataFrame containing KW test results.
  70. df_pval_corrected : pandas.DataFrame
  71. DataFrame containing Dunn's test results with corrected p-values.
  72. significant_syllables : list
  73. List of corrected KW significant syllables (syllables with p-values < thresh).
  74. """
  75. rnd = np.random.RandomState(seed=seed)
  76. # get grouped mean data
  77. grouped_data = (
  78. stats_df.pivot_table(
  79. index=["group", "name"], columns="syllable", values=statistic
  80. )
  81. .replace(np.nan, 0)
  82. .reset_index()
  83. )
  84. # compute KW constants
  85. vc = grouped_data.group.value_counts().loc[grouped_data.group.unique()]
  86. n_per_group = vc.values
  87. group_names = vc.index
  89. cum_group_idx = np.insert(np.cumsum(n_per_group), 0, 0)
  90. num_groups = len(group_names)
  92. # get all syllable usage data
  93. df_only_stats = grouped_data.drop(["group", "name"], axis=1)
  94. syllable_data = grouped_data.drop(["group", "name"], axis=1).values
  96. N_m, N_s = syllable_data.shape
  97. # Run KW and return H-stats
  98. h_all, real_ranks, X_ties = run_manual_KW_test(
  99. df_usage=df_only_stats,
  100. merged_usages_all=syllable_data,
  101. num_groups=num_groups,
  102. n_per_group=n_per_group,
  103. cum_group_idx=cum_group_idx,
  104. n_perm=n_perm,
  105. seed=seed,
  106. )
  108. # find the real k_real
  109. df_k_real = pd.DataFrame(
  110. [
  111. stats.kruskal(
  112. *np.array_split(syllable_data[:, s_i], np.cumsum(n_per_group[:-1]))
  113. )
  114. for s_i in range(N_s)
  115. ]
  116. )
  118. # multiple test correction
  119. df_k_real["p_adj"] = multipletests(
  120. ((h_all > df_k_real.statistic.values).sum(0) + 1) / n_perm,
  121. alpha=thresh,
  122. method=mc_method,
  123. )[1]
  125. # return significant syllables based on the threshold
  126. df_k_real["is_sig"] = df_k_real["p_adj"] <= thresh
  128. # Run Dunn's z-test statistics
  129. (
  130. null_zs_within_group,
  131. real_zs_within_group,
  132. ) = dunns_z_test_permute_within_group_pairs(
  133. grouped_data, vc, real_ranks, X_ties, N_m, group_names, rnd, n_perm
  134. )
  136. # Compute p-values from Dunn's z-score statistics
  137. df_pair_corrected_pvalues, _ = compute_pvalues_for_group_pairs(
  138. real_zs_within_group,
  139. null_zs_within_group,
  140. df_k_real,
  141. group_names,
  142. n_perm,
  143. thresh,
  144. mc_method,
  145. )
  147. # combine Dunn's test results into single DataFrame
  148. df_z = pd.DataFrame(real_zs_within_group)
  149. df_z.index = df_z.index.set_names("syllable")
  150. dunn_results_df = df_z.reset_index().melt(id_vars=[("syllable", "")])
  151. dunn_results_df.rename(
  152. columns={"variable_0": "group1", "variable_1": "group2"}, inplace=True
  153. )
  155. # Get intersecting significant syllables between
  156. intersect_sig_syllables = {}
  157. pvals = {}
  158. for pair in df_pair_corrected_pvalues.columns.tolist():
  159. intersect_sig_syllables[pair] = np.where(
  160. (df_pair_corrected_pvalues[pair] < thresh) & (df_k_real.is_sig)
  161. )[0]
  162. pvals[pair] = df_pair_corrected_pvalues[pair]
  164. return df_k_real, dunn_results_df, intersect_sig_syllables, pvals
  166. def run_manual_KW_test(
  167. df_usage,
  168. merged_usages_all,
  169. num_groups,
  170. n_per_group,
  171. cum_group_idx,
  172. n_perm=10000,
  173. seed=42,
  174. ):
  175. """Run a manual Kruskal-Wallis test compare the results agree with the
  176. scipy.stats.kruskal function.
  178. Parameters
  179. ----------
  180. df_usage : pandas.DataFrame
  181. DataFrame with syllable usages. shape = (N_m, n_syllables)
  182. merged_usages_all : np.array
  183. numpy array format of the df_usage DataFrame.
  184. num_groups : int
  185. Number of unique groups
  186. n_per_group : list
  187. list of value counts for recordings per group. len == num_groups.
  188. cum_group_idx : list
  189. list of indices for different groups. len == num_groups + 1.
  190. n_perm : int, optional
  191. Number of permuted samples to generate, by default 10000
  192. seed : int, optional
  193. Random seed used to initialize the pseudo-random number generator, by default 42
  195. Returns
  196. -------
  197. h_all : np.array
  198. Array of H-stats computed for given n_syllables; shape = (n_perms, N_s)
  199. real_ranks : np.array
  200. Array of syllable ranks, shape = (N_m, n_syllables)
  201. X_ties : np.array
  202. 1-D list of tied ranks, where if value > 0, then rank is tied. len(X_ties) = n_syllables
  203. """
  205. N_m, N_s = merged_usages_all.shape
  207. # create random index array n_perm times
  208. rnd = np.random.RandomState(seed=seed)
  209. perm = rnd.rand(n_perm, N_m).argsort(-1)
  211. # get degrees of freedom
  212. dof = num_groups - 1
  214. real_ranks = np.apply_along_axis(stats.rankdata, 0, merged_usages_all)
  215. X_ties = df_usage.apply(get_tie_correction, 0, N_m=N_m).values
  216. KW_tie_correct = np.apply_along_axis(stats.tiecorrect, 0, real_ranks)
  218. # rank data
  219. perm_ranks = real_ranks[perm]
  221. # get square of sums for each group
  222. ssbn = np.zeros((n_perm, N_s))
  223. for i in range(num_groups):
  224. ssbn += (
  225. perm_ranks[:, cum_group_idx[i] : cum_group_idx[i + 1]].sum(1) ** 2
  226. / n_per_group[i]
  227. )
  229. # h-statistic
  230. h_all = 12.0 / (N_m * (N_m + 1)) * ssbn - 3 * (N_m + 1)
  231. h_all /= KW_tie_correct
  232. p_vals = stats.chi2.sf(h_all, df=dof)
  234. # check that results agree
  235. p_i = np.random.randint(n_perm)
  236. s_i = np.random.randint(N_s)
  237. kr = stats.kruskal(
  238. *np.array_split(
  239. merged_usages_all[perm[p_i, :], s_i], np.cumsum(n_per_group[:-1])
  240. )
  241. )
  242. assert (kr.statistic == h_all[p_i, s_i]) & (
  243. kr.pvalue == p_vals[p_i, s_i]
  244. ), "manual KW is incorrect"
  246. return h_all, real_ranks, X_ties
  248. def get_tie_correction(x, N_m):
  249. """Assign tied rank values to the average of the ranks they would have
  250. received if they had not been tied for Kruskal-Wallis helper function.
  252. Parameters
  253. ----------
  254. x : pd.Series
  255. syllable usages for a single recording.
  256. N_m : int
  257. Number of total recordings.
  259. Returns
  260. -------
  261. corrected_rank : float
  262. average of the inputted tied ranks.
  263. """
  265. vc = x.value_counts()
  266. tie_sum = 0
  267. if (vc > 1).any():
  268. tie_sum += np.sum(vc[vc != 1] ** 3 - vc[vc != 1])
  269. return tie_sum / (12.0 * (N_m - 1))
  271. def dunns_z_test_permute_within_group_pairs(
  272. df_usage, vc, real_ranks, X_ties, N_m, group_names, rnd, n_perm
  273. ):
  274. """Run Dunn's z-test statistic on combinations of all group pairs, handling
  275. pre- computed tied ranks.
  277. Parameters
  278. ----------
  279. df_usage : pandas.DataFrame
  280. DataFrame containing only pre-computed syllable stats.
  281. vc : pd.Series
  282. value counts of recordings in each group.
  283. real_ranks : np.array
  284. Array of syllable ranks.
  285. X_ties : np.array
  286. 1-D list of tied ranks, where if value > 0, then rank is tied
  287. N_m : int
  288. Number of recordings.
  289. group_names : pd.Index
  290. Index list of unique group names.
  291. rnd : np.random.RandomState
  292. Pseudo-random number generator.
  293. n_perm : int
  294. Number of permuted samples to generate.
  296. Returns
  297. -------
  298. null_zs_within_group : dict
  299. dict of group pair keys paired with vector of Dunn's z-test statistics of the null hypothesis.
  300. real_zs_within_group : dict
  301. dict of group pair keys paired with vector of Dunn's z-test statistics
  302. """
  304. null_zs_within_group = {}
  305. real_zs_within_group = {}
  307. A = N_m * (N_m + 1.0) / 12.0
  309. for i_n, j_n in combinations(group_names, 2):
  310. is_i = df_usage.group == i_n
  311. is_j = df_usage.group == j_n
  313. n_mice = is_i.sum() + is_j.sum()
  315. ranks_perm = real_ranks[(is_i | is_j)][rnd.rand(n_perm, n_mice).argsort(-1)]
  316. diff = np.abs(
  317. ranks_perm[:, : is_i.sum(), :].mean(1)
  318. - ranks_perm[:, is_i.sum() :, :].mean(1)
  319. )
  320. B = 1.0 / vc.loc[i_n] + 1.0 / vc.loc[j_n]
  322. # also do for real data
  323. group_ranks = real_ranks[(is_i | is_j)]
  324. real_diff = np.abs(
  325. group_ranks[: is_i.sum(), :].mean(0) - group_ranks[is_i.sum() :, :].mean(0)
  326. )
  328. # add to dict
  329. pair = (i_n, j_n)
  330. null_zs_within_group[pair] = diff / np.sqrt((A - X_ties) * B)
  331. real_zs_within_group[pair] = real_diff / np.sqrt((A - X_ties) * B)
  333. return null_zs_within_group, real_zs_within_group
  335. def compute_pvalues_for_group_pairs(
  336. real_zs_within_group,
  337. null_zs,
  338. df_k_real,
  339. group_names,
  340. n_perm=10000,
  341. thresh=0.05,
  342. mc_method="fdr_bh",
  343. ):
  344. """Adjust the p-values from Dunn's z-test statistics and computes the
  345. resulting significant syllables with the adjusted p-values.
  347. Parameters
  348. ----------
  349. real_zs_within_group : dict
  350. dict of group pair keys paired with vector of Dunn's z-test statistics
  351. null_zs : dict
  352. dict of group pair keys paired with vector of Dunn's z-test statistics of the null hypothesis.
  353. df_k_real : pandas.DataFrame
  354. DataFrame of KW test results.
  355. group_names : pd.Index
  356. Index list of unique group names.
  357. n_perm : int, optional
  358. Number of permuted samples to generate, by default 10000
  359. thresh : float, optional
  360. Alpha threshold to consider syllable significant, by default 0.05
  361. mc_method : str, optional
  362. Multiple Corrections method to use, by default "fdr_bh"
  363. verbose : bool, optional
  364. indicates whether to print out the significant syllable results, by default False
  366. Returns
  367. -------
  368. df_pval_corrected : pandas.DataFrame
  369. DataFrame containing Dunn's test results with corrected p-values.
  370. significant_syllables : list
  371. List of corrected KW significant syllables (syllables with p-values < thresh).
  372. """
  374. # do empirical p-val calculation for all group permutation
  376. p_vals_allperm = {}
  377. for pair in combinations(group_names, 2):
  378. p_vals_allperm[pair] = (
  379. (null_zs[pair] > real_zs_within_group[pair]).sum(0) + 1
  380. ) / n_perm
  382. # summarize into df
  383. df_pval = pd.DataFrame(p_vals_allperm)
  385. def correct_p(x):
  386. return multipletests(x, alpha=thresh, method=mc_method)[1]
  388. df_pval_corrected = df_pval.apply(correct_p, axis=1, result_type="broadcast")
  390. return df_pval_corrected, ((df_pval_corrected[df_k_real.is_sig] < thresh).sum(0))
  392. def sort_syllables_by_stat(stats_df, stat="frequency"):
  393. """Sort sylllabes by the stat and return the ordering and label mapping.
  395. Parameters
  396. ----------
  397. stats_df : pandas.DataFrame
  398. the stats dataframe that contains kinematic data and the syllable label for each recording and each syllable
  399. stat : str, optional
  400. the statistic to sort on, by default 'frequency'
  402. Returns
  403. -------
  404. ordering : list
  405. the list of syllables sorted by the stat
  406. relabel_mapping : dict
  407. the mapping from the syllable to the new plotting label
  408. """
  410. # stats_df frequency normalized by session
  411. # mean frequency by syllable don't always refect the ordering from reindexing
  412. # use the syllable label as ordering instead
  413. if stat == "frequency":
  414. ordering = sorted(stats_df.syllable.unique())
  415. else:
  416. ordering = (
  417. stats_df.drop(
  418. [col for col, dtype in stats_df.dtypes.items() if dtype == "object"],
  419. axis=1,
  420. )
  421. .groupby("syllable")
  422. .mean()
  423. .sort_values(by=stat, ascending=False)
  424. .index
  425. )
  427. # Get sorted ordering
  428. ordering = list(ordering)
  430. # Get order mapping
  431. relabel_mapping = {o: i for i, o in enumerate(ordering)}
  433. return ordering, relabel_mapping
  435. # %%
  436. # read in sig_syllables if previously saved into csv
  438. # sig_syllables = []
  439. # all_sig_syllables = pd.read_csv(project_dir+model_name+'/'+sig_syllables_filename)
  440. # for i in range(all_sig_syllables.shape[0]):
  441. # if all_sig_syllables.iloc[i,0] == group_1 and all_sig_syllables.iloc[i,1] == group_2:
  442. # sig_syllables = all_sig_syllables.iloc[i,2]
  443. # elif all_sig_syllables.iloc[i,0] == group_2 and all_sig_syllables.iloc[i,1] == group_1:
  444. # sig_syllables = all_sig_syllables.iloc[i,2]
  445. # sig_syllables
  447. # %%
  448. # extract significant pairs and p-values from stats_df
  449. _, _, sig_pairs, pvals = run_kruskal(stats_df)
  451. # %%
  452. # save pairwise p-vals in csv
  453. import csv
  455. with open(os.path.join(stats_dir,'pairwise_pvals.csv'), 'w') as csvfile:
  456. writer = csv.writer(csvfile)
  458. for item in pvals:
  459. first_group = item[0]
  460. second_group = item[1]
  461. pval = (pvals[item])
  462. writer.writerow([first_group,second_group,pval])
  464. # %%
  465. # specify the two groups for pairwise comparison and plot relative frequencies of syllables, marking significant syllables with a red star
  466. # right over the x axis
  468. current_stats_df = stats_df[(stats_df['group'] == group_1) | (stats_df['group'] == group_2)]
  470. ordering, _ = sort_syllables_by_stat(stats_df, stat=stat)
  472. if (group_1, group_2) in sig_pairs.keys():
  473. sig_sylls = sig_pairs.get((group_1, group_2))
  474. else:
  475. sig_sylls = sig_pairs.get((group_2, group_1))
  477. hue = 'group'
  479. fig, ax = plt.subplots(1,1,figsize=(6,4))
  480. ax = sns.pointplot(
  481. data=current_stats_df,
  482. x="syllable",
  483. y=stat,
  484. hue=hue,
  485. palette='Reds',
  486. linestyle="none",
  487. errorbar=('se'),
  488. markersize=3,
  489. marker="o",
  490. err_kws={'linewidth':1}
  491. )
  492. syllables = stats_df["syllable"].unique()
  493. markings = []
  494. for s in sig_sylls:
  495. if s in ordering:
  496. markings.append(np.where(ordering == s)[0])
  497. else:
  498. continue
  499. if len(markings) > 0:
  500. markings = np.concatenate(markings)
  501. plt.scatter(markings, [0] * len(markings), color="r", marker="*",linewidths=0.05)
  502. else:
  503. print("No significant syllables found.")
  505. ax.spines["top"].set_visible(False)
  506. ax.spines["right"].set_visible(False)
  507. ax.spines["left"].set_linewidth(0.5)
  508. ax.spines["bottom"].set_linewidth(0.5)
  509. ax.legend(frameon=False, loc='upper right')
  510. ax.set_ylim([-0.003,None])
  512. plt.tight_layout()
  514. # save fig as png and eps in figs directory
  515. fig.savefig(project_dir+model_name+'/'+fig_dir+'p_'+group_1+'_'+group_2+'_syllable_freq.png')
  516. fig.savefig(project_dir+model_name+'/'+fig_dir+'p_'+group_1+'_'+group_2+'_syllable_freq.eps', format='eps')
  519. # %%
  520. # plot correlation graphs for a phenotype of interest, here avoidance index (ai)
  521. # ai data was saved in an xlsx file with the columns 'file' (subject) and 'ai'
  523. for s in range(len(stats_df["syllable"].unique())):
  524. current_syllable = s
  525. current_stats_df = stats_df[(stats_df['syllable'] == current_syllable)]
  526. correlation_df = pd.DataFrame(columns=['group','frequency','ai'])
  527. # insert name of excel file with phenotype data (2 columns, indicate column names in usecols=[])
  528. ai_df = pd.read_excel('',usecols=['file','ai'])
  529. for i in range(current_stats_df.shape[0]):
  530. new_entry_df = pd.DataFrame(columns=['group','frequency','ai'],index=range(1))
  531. group = current_stats_df.iloc[i,0]
  532. name = current_stats_df.iloc[i,1]
  533. frequency = current_stats_df.iloc[i,15]
  534. # find ai for current subject in ai_df
  535. for j in range(ai_df.shape[0]):
  536. if name == ai_df.iloc[j,0]:
  537. ai = ai_df.iloc[j,1]
  538. new_entry_df.iloc[0,0] = group
  539. new_entry_df.iloc[0,1] = float(frequency)
  540. new_entry_df.iloc[0,2] = float(ai)
  541. correlation_df = pd.concat([correlation_df,new_entry_df])
  543. correlation_df = correlation_df.astype({'frequency':'float','ai':'float'})
  545. color = sns.color_palette(palette='BuPu')[5]
  547. # creates correlation plots colored by group
  548. lm = sns.lmplot(data=correlation_df,x='frequency',y='ai',col="group",facet_kws=dict(sharex=False, sharey=False),scatter_kws={"color":color}, line_kws={"color":color})
  549. font={'fontname':'Arial'}
  550. for group, ax in lm.axes_dict.items():
  551. r, p = stats.pearsonr(correlation_df[correlation_df['group'] == group]['frequency'],correlation_df[correlation_df['group'] == group]['ai'])
  552. r2=r**2
  553. ax.text(0.05,1,'$r^2$={:.2f}, p={:.2g}'.format(r2,p),horizontalalignment='left',verticalalignment='top',transform=ax.transAxes)
  554. ax.spines["top"].set_visible(False)
  555. ax.spines["right"].set_visible(False)
  556. ax.spines["left"].set_linewidth(0.5)
  557. ax.spines["bottom"].set_linewidth(0.5)
  558. ax.set_ylim([-1,1])
  559. plt.suptitle('syllable '+str(current_syllable),fontsize=20,y=1.05)
  561. lm.savefig(project_dir+model_name+'/'+fig_dir+'syllable '+str(current_syllable)+'_corr_by_group.png')
  563. for group, ax in lm.axes_dict.items():
  564. plt.clf()
  565. plt.close()
  567. # creates correlation plots of all combined data
  568. rplot = sns.regplot(data=correlation_df,x='frequency',y='ai',scatter_kws={"color":color}, line_kws={"color":color})
  569. ax = rplot.axes
  570. r, p = stats.pearsonr(correlation_df['frequency'],correlation_df['ai'])
  571. r2=r**2
  572. ax.text(0.05,1,'$r^2$={:.2f}, p={:.2g}'.format(r2,p),horizontalalignment='left',verticalalignment='top',transform=ax.transAxes)
  573. ax.spines["top"].set_visible(False)
  574. ax.spines["right"].set_visible(False)
  575. ax.spines["left"].set_linewidth(0.5)
  576. ax.spines["bottom"].set_linewidth(0.5)
  577. ax.set_ylim([-1,1])
  578. plt.title('syllable '+str(current_syllable),fontsize=16,y=1.05)
  579. plt.tight_layout()
  580. figure = rplot.get_figure()
  581. figure.savefig(project_dir+model_name+'/'+fig_dir+'syllable '+str(current_syllable)+'_corr_combined.png')
  582. plt.clf()
  583. plt.close()
  585. # %%
  586. sns.color_palette(palette='BuPu')
  588. # %%
  589. # format of moseq_df
  590. moseq_df
  592. # %%
  593. # for analysis with respect to the subject's position in the video, can create a csv indicating boundaries
  594. # here, x positions are defined for each video to determine whether the mouse is in the left chamber or right chamber
  595. # if arena does not move between videos, can be manually entered
  596. # here, since arena moves between videos, x cutoffs are defined for each video/subject
  598. # read in left and right x cutoff values for each video
  599. l_r_df = pd.read_csv('/Users/claraliff/Desktop/Lab/eLife_paper/kpMoSeq/v3/v3_f0_f1_moseq/2024_11_18-12_14_41/l_r_xcutoffs.csv')
  600. # define whether left or right is conditioned odor side using an experiment metadata excel sheet
  601. c_o_side = pd.read_excel('/Users/claraliff/Desktop/Lab/eLife_paper/kpMoSeq/v3/v3_f0_f1_csvs_grouped.xlsx',usecols="A,D")
  602. port_x_df = pd.DataFrame(columns=['name','port_x'],index=range(99))
  603. port_x_df.loc[:,'name'] = c_o_side['file']
  604. port_x_df
  606. # in our videos, the conditioned odor port, a region of interest in our videos, was always located 200 pixels to the left of the left
  607. # chamber entrance or 200 pixels to the right of the right chamber entrance
  608. # port_x_df defines the x position of the conditioned odor port for each video/subject, to later be used to calculate real-time distance of
  609. # the mouse from the odor source
  611. for x in range(l_r_df.shape[0]):
  612. if c_o_side.iloc[x,1] == "R":
  613. port_x = l_r_df.iloc[x,2] + 200
  614. elif c_o_side.iloc[x,1] == "L":
  615. port_x = l_r_df.iloc[x,1] - 200
  616. port_x_df.iloc[x,1] = port_x
  618. port_x_df
  620. # %%
  621. # append two columns to moseq_df: the x position of the conditioned odor port, and the x distance of the mouse from the odor port
  623. moseq_df_dist = moseq_df.copy()
  624. moseq_df_dist["port_x"] = np.nan
  625. moseq_df_dist["port_dist"] = np.nan
  626. moseq_df_dist
  628. # %%
  629. # for each row (frame) of all videos, calculate the x distance of the mouse x centroid from the x coordinate of the conditioned odor port
  631. for x in range(moseq_df_dist.shape[0]):
  632. current_name = moseq_df_dist.iloc[x,0]
  633. i = port_x_df.index[port_x_df.iloc[:,0].str.contains(current_name)]
  634. port_x = port_x_df.iloc[i,1]
  635. moseq_df_dist.iloc[x,10] = port_x
  636. current_centroid_x = moseq_df_dist.iloc[x,1]
  637. port_dist = abs(port_x - current_centroid_x)
  638. moseq_df_dist.iloc[x,11] = port_dist
  639. moseq_df_dist
  641. # %%
  642. # sanity check that all distances are positive
  643. moseq_df_dist[moseq_df_dist['port_dist']<0]
  645. # %%
  646. # save moseq_df with distance info
  647. moseq_df_dist.to_csv(project_dir+model_name+'/moseq_df/moseq_df_dist.csv')
  649. # %%
  650. # read in moseq_df_dist if continuing after checkpoint
  651. # moseq_df_dist = pd.read_csv(project_dir+model_name+'/moseq_df/moseq_df_dist.csv')
  653. # %%
  654. # create temp_df, a subset of moseq_df_dist with your two groups for comparison
  655. # filtering to only analyze top 20 (20 most frequent) syllables
  657. temp_df = moseq_df_dist[(moseq_df_dist['group'] == 'f1_up') | (moseq_df_dist['group'] == 'f1_p')]
  658. temp_df = temp_df[temp_df['syllable'] < 20]
  659. temp_df
  661. # %%
  662. # create variables for x position (all frames, all subjects) for each group
  663. f1_p_x_dist = temp_df[temp_df['group'] == 'f1_p']
  664. f1_p_x_dist = f1_p_x_dist.iloc[:,2]
  665. f1_p_x_dist
  667. f1_up_x_dist = temp_df[temp_df['group'] == 'f1_up']
  668. f1_up_x_dist = f1_up_x_dist.iloc[:,2]
  669. f1_up_x_dist
  671. # %%
  672. # plot kde for x position by group
  673. # black bars indicate the left/right chamber entries (cutoffs)
  675. d = sns.displot(temp_df,x='port_dist',palette='BuPu',hue='group',kind="kde",common_norm=True)
  676. d.set(xlim=(0,600))
  677. plt.axvline(200,0,1,color='black')
  678. plt.axvline(350,0,1,color='black')
  679. d.savefig(project_dir+model_name+'/'+fig_dir+'port_dist_hist_kde_common_norm.eps', format='eps')
  681. # %%
  682. # plot as histogram
  684. d = sns.displot(temp_df,x='port_dist',palette='BuPu',hue='group',kind="hist",stat='frequency',bins=50)
  685. d.set(xlim=(0,600))
  686. plt.axvline(200,0,1,color='black')
  687. plt.axvline(350,0,1,color='black')
  688. # d.savefig(project_dir+model_name+'/'+fig_dir+'port_dist_hist_freq_bins25.eps', format='eps')
  689. # d.savefig(project_dir+model_name+'/'+fig_dir+'port_dist_hist_freq_bins25.png', format='png')
  691. # %%
  692. # plot syllable usage across relative x position for each syllable
  693. dp = sns.displot(temp_df,x='port_dist',palette='BuPu',col='syllable',hue='group',col_wrap=5,facet_kws={'sharey':False},kind='kde',common_norm=False)
  694. for ax in dp.axes.flat:
  695. ymin,ymax=ax.get_ylim()
  696. ax.set(xlim=(0,600))
  697. ax.vlines(200,ymin,ymax,colors='black')
  698. ax.vlines(350,ymin,ymax,colors='black')
  699. dp.savefig(project_dir+model_name+'/'+fig_dir+'port_dist_hist_by_syll.eps', format='eps')

plot_kpmoseq_data.ipynb at commit cc8935b, under MIT · at the source

Overview

Authors: Clara W Liff1,2, Yasmine R Ayman1, Eliza CB Jaeger1, Avery Cardeiro1, Hudson S Lee1, Alexis Kim1, Angelica Vina-Abarracin1, Dianne-Lee KD Ferguson1,3, Bianca J Marlin1,2,3,4
ORCID iDs: Clara W Liff, Bianca J Marlin
  1. Mortimer B. Zuckerman Mind Brain and Behavior Institute, Columbia University, New York, United States
  2. Department of Neuroscience, Columbia University, New York, United States
  3. Howard Hughes Medical Institute, Columbia University, New York, United States
  4. Department of Psychology, Columbia University, New York, United States
Institutions: Mortimer B. Zuckerman Mind Brain Behavior Institute (United States); Columbia University (United States); Howard Hughes Medical Institute (United States)
Journal: eLife, volume 12, article RP92882
Dates: published online 14 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.92882 · PMID 41979323 · PMCID PMC13078775 · OpenAlex W4389781873
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), cellular / molecular (subfield)
Methods: Statistics, Evoked potentials, Machine learning
Keywords: Mouse
MeSH: Conditioning, Psychological*, Fear*, Olfactory Mucosa*, Olfactory Receptor Neurons*, Animals, Female, Male, Mice, Mice, Inbred C57BL, Odorants, Smell (* major topic)
Topic: Olfactory and Sensory Function Studies (Sensory Systems, Neuroscience), according to OpenAlex
Funding: United Negro College Fund (E.E. Just Fellowship CU20-1071); Howard Hughes Medical Institute (Freeman Hrabrowski Scholar); NIH HHS (1S10OD023587-01); Simons Foundation (Simons Society of Fellows Junior Fellowship 524991); NIMH NIH HHS (T32MH126036); Brain and Behavior Research Foundation (NARSAD Young Investigator Grant 30380)
Citations: cited by 3 papers (Europe PMC); 50 references in the paper
Research resources: RRID:AB_10000240, RRID:AB_2535866

Abstract

The main olfactory epithelium initiates the process of odor encoding. Recent studies have demonstrated intergenerationally inherited changes in the olfactory system in response to fear conditioning, resulting in increases in olfactory sensory neuron frequencies and altered responses to odors. We investigated changes in the cellular composition of the olfactory epithelium in response to an aversive stimulus. Here, we achieve volumetric cellular resolution to demonstrate that olfactory fear conditioning increases the number of odor-encoding neurons in mice that experience odor-shock conditioning (F0), as well as their unconditioned offspring (F1). We demonstrate that the increase in F0 is due, in part, to the biasing of the stem cell layer of the main olfactory epithelium. A detailed analysis of F1 behavior revealed subtle odor-specific differences between the offspring of unconditioned and conditioned parents, despite the absence of an active aversion to the conditioned odor. Thus, we reveal intergenerational regulation of olfactory epithelium composition in response to olfactory fear conditioning, providing insight into the heritability of acquired phenotypes.

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

Repository

Its files are read in the Code ↔ Paper reader above.

BJMarlinLab/Liff_et_al_2026

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: cc8935bcac1028980acff49a03a51e546746420b, 25 March 2026
Languages: Jupyter (1)
Size: 3 files, 1 script
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 1 notebook
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
3 files
At the source: github.com/BJMarlinLab/Liff_et_al_2026

Code availability

Code is available at https://github.com/BJMarlinLab/Liff_et_al_2026 (copy archived at Marlin Lab, 2026 (https://github.com/BJMarlinLab/Liff_et_al_2026)).

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

Tracing map

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

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • no match between paragraphs and code yet;
  • 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

Dataset available on Dryad at DOI: https://doi.org/10.5061/dryad.80gb5mm4m. The M71-RFP mouse line is available upon request.

The following dataset was generated:

Liff C, Ayman Y, Jaeger E, Cardeiro A, Lee H, Kim A, Vina-Albarracin A, Ferguson D-L, Marlin B. 2026. Fear conditioning biases olfactory sensory neuron frequencies across generations. Dryad Digital Repository.

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

Recorded: type, language, journal, volume, pages, dates, 9 authors, 1 keyword, 11 MeSH terms, 6 funders, 49 references, 2 RRIDs.

Cite

This paper

Liff, C. W., Ayman, Y. R., Jaeger, E. C., Cardeiro, A., Lee, H. S., Kim, A., Vina-Abarracin, A., Ferguson, D.-L. K., & Marlin, B. J. (2026). Fear conditioning biases olfactory sensory neuron frequencies across generations. eLife, 12, RP92882. https://doi.org/10.7554/elife.92882

BibTeX

@article{liff2026fear,
author = {Liff, Clara W and Ayman, Yasmine R and Jaeger, Eliza CB and Cardeiro, Avery and Lee, Hudson S and Kim, Alexis and Vina-Abarracin, Angelica and Ferguson, Dianne-Lee KD and Marlin, Bianca J},
title = {{Fear conditioning biases olfactory sensory neuron frequencies across generations}},
journal = {eLife},
year = {2026},
month = apr,
volume = {12},
pages = {RP92882},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.92882},
url = {https://doi.org/10.7554/elife.92882},
pmid = {41979323},
pmcid = {PMC13078775}
}

RIS

TY - JOUR
AU - Liff, Clara W
AU - Ayman, Yasmine R
AU - Jaeger, Eliza CB
AU - Cardeiro, Avery
AU - Lee, Hudson S
AU - Kim, Alexis
AU - Vina-Abarracin, Angelica
AU - Ferguson, Dianne-Lee KD
AU - Marlin, Bianca J
TI - Fear conditioning biases olfactory sensory neuron frequencies across generations
T2 - eLife
J2 - Elife
PY - 2026
DA - 2026/04/14
VL - 12
SP - RP92882
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.92882
UR - https://doi.org/10.7554/elife.92882
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.92882",
"type": "article-journal",
"title": "Fear conditioning biases olfactory sensory neuron frequencies across generations",
"container-title": "eLife",
"author": [
{
"family": "Liff",
"given": "Clara W"
},
{
"family": "Ayman",
"given": "Yasmine R"
},
{
"family": "Jaeger",
"given": "Eliza CB"
},
{
"family": "Cardeiro",
"given": "Avery"
},
{
"family": "Lee",
"given": "Hudson S"
},
{
"family": "Kim",
"given": "Alexis"
},
{
"family": "Vina-Abarracin",
"given": "Angelica"
},
{
"family": "Ferguson",
"given": "Dianne-Lee KD"
},
{
"family": "Marlin",
"given": "Bianca J"
}
],
"container-title-short": "Elife",
"volume": "12",
"page": "RP92882",
"DOI": "10.7554/elife.92882",
"PMID": "41979323",
"PMCID": "PMC13078775",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.92882",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
14
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s41467-026-72445-1 [code]
Perception and neural representation of intermittent odor stimuli in mice.
Journal: Nature communications
In common: statsmodels, seaborn, pandas, 3 other tools, mouse, 3 references
[2] doi:10.1038/s41598-026-51082-0
Non-uniform endogenous regeneration of olfactory sensory neuron axons across the mouse olfactory bulb.
Journal: Scientific reports
In common: mouse, 5 references
[3] doi:10.1371/journal.pgen.1012090 [code]
Delta family protocadherins contribute to protoglomerular targeting of olfactory sensory neuron axons in the olfactory bulb.
Journal: PLoS genetics
In common: cellular / molecular, 5 references
[4] doi:10.1038/s41467-026-71595-6 [code]
A single-cell and spatial atlas of early human olfactory development.
Journal: Nature communications
In common: seaborn, pandas, Matplotlib, 1 other tool, 3 references
[5] doi:10.1038/s41586-026-10444-4 [code]
A brain reward circuit inhibited by next-generation weight-loss drugs in mice.
Journal: Nature
In common: seaborn, pandas, SciPy, 2 other tools, mouse, 2 references
[6] doi:10.1038/s41593-026-02262-8 [code]
Cheese3D enables sensitive detection and analysis of whole-face movement in mice.
Journal: Nature neuroscience
In common: seaborn, pandas, SciPy, 2 other tools, mouse, 2 references
[7] doi:10.1016/j.celrep.2026.117419 [code]
Conserved role of primary motor cortex in the control of prehension in mice and macaques.
Journal: Cell reports
In common: statsmodels, seaborn, pandas, 3 other tools, mouse, 1 reference
[8] 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: statsmodels, seaborn, pandas, 3 other tools, 1 reference
[9] doi:10.1038/s41586-026-10679-1 [code]
Cortical development dynamics across autism spectrum disorder mouse models.
Journal: Nature
In common: seaborn, pandas, SciPy, 2 other tools, mouse, cellular / molecular, 1 reference
[10] doi:10.1016/j.cub.2026.06.016 [code]
Neuronal RNAi and oxygen-sensing circuit shape germline resilience to heat stress.
Journal: Current biology : CB
In common: statsmodels, pandas, SciPy, 2 other tools, 1 reference

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.