OSCR

The gut microbiota and sleep in infants: a focus on diurnal rhythmicity patterns.

Code ↔ Paper

16 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 16 matches · 2 of them tie a paragraph to a whole file without given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Data processing › Infant gut microbiota › Sequences processing: gut microbiota composition and diversity. ↔ scripts/functions_script.py, lines 16 to 43 · score 0.86 · Faith phylogenetic diversity, Bray Curtis dissimilarity, Shannon entropy, unweighted UniFrac, Pielou evenness, Jaccard · Report a mismatch
  2. [2] § Methods › Data processing › Infant sleep › Sleep quality: infant sleep quality index (BabySQUID). ↔ scripts/2-EDA.ipynb, lines 137 to 154 · score 0.82 · sleep onset latencies, nighttime sleep duration, nighttime awakenings, babySQUID, BISQ, bedtime · Report a mismatch
  3. [3] § Results › Modulation of infant sleep by the gut microbiota › Sleep quality: BabySQUID ↔ scripts/2-EDA.ipynb, lines 137 to 154 · score 0.81 · sleep onset latency, nighttime sleep duration, nighttime awakenings, infant sleep BISQ, babySQUID, bedtime · Report a mismatch
  4. [4] § Methods › Data processing › Infant sleep › Sleep quality: infant sleep quality index (BabySQUID). ↔ scripts/functions_script.py, lines 347 to 385 · score 0.81 · min max, nighttime sleep duration, equal weights, reverse, awakenings, latency · Report a mismatch
  5. [5] § Methods › Statistical analysis › Modulation of infant sleep by the gut microbiota › Sleep quality: BabySQUID. ↔ scripts/functions_script.py, lines 1336 to 1448 · score 0.76 · random forest classifiers, cross validation, Accuracy, splits, trained, predict · Report a mismatch
  6. [6] § Methods › Data processing › Infant gut microbiota › Sequences processing: gut microbiota composition and diversity. ↔ scripts/1-seq-proc.ipynb, lines 116 to 162 · score 0.76 · classifier trained, sequencing depth, bp, RESCRIPt, SILVA, truncation · Report a mismatch
  7. [7] § Results › Modulation of infant sleep by the gut microbiota › Sleep rhythmicity: CFI ↔ scripts/functions_script.py, lines 16 to 43 · score 0.66 · Faith phylogenetic diversity, attuned care, Bray Curtis dissimilarity, temporal volatility, Jaccard, alpha diversity · Report a mismatch
  8. [8] § Methods › Statistical analysis › Exploratory analysis of the samples › Gut melatonin. ↔ scripts/6-melatonin.ipynb, lines 56 to 65 · score 0.65 · bowel movement, spent awake, gut melatonin, pg, stool · Report a mismatch
  9. [9] § Methods › Statistical analysis › Modulation of infant sleep by the gut microbiota › Sleep quality: BabySQUID. ↔ scripts/5.3-mb-sleep-RF.ipynb, lines 133 to 157 · score 0.64 · cross validation, splits, ASV, trained, harmonic, precision · Report a mismatch
  10. [10] § Results › Diurnal patterns along age › Gut microbiota rhythmicity and sleep rhythmicity ↔ scripts/1.2-meta-proc.ipynb, lines 402 to 426 · score 0.59 · Clostridium sensu stricto, Escherichia Shigella, Bacteroides · Report a mismatch
  11. [11] § Methods › Statistical analysis › Exploratory analysis of the samples › Gut melatonin. ↔ scripts/qiime2/14-diff-abund.sh, the whole file · a weak match · score 0.57 · ANCOM BC2, prevalence, filter, QIIME, abundance, melatonin · Report a mismatch
  12. [12] § Results › Modulation of infant sleep by the gut microbiota › Sleep rhythmicity: CFI ↔ scripts/functions_script.py, lines 45 to 57 · score 0.57 · Shannon entropy, Pielou evenness, phylogenetic diversity, temporal volatility, alpha diversity, Faith · Report a mismatch
  13. [13] § Results › Diurnal patterns along age › Microbial temporal volatility and age ↔ scripts/4-factors.ipynb, lines 93 to 116 · score 0.55 · Bray Curtis dissimilarity, unweighted UniFrac, Temporal volatility, distance, age · Report a mismatch
  14. [14] § Methods › Data processing › Infant gut microbiota › Stool processing: DNA extraction and 16S rRNA gene sequencing. ↔ scripts/qiime2/4-get_silva_reads_n_classifier.sh, the whole file · a weak match · score 0.54 · DNA, primers, V4, reverse, Amplicon, sequencing · Report a mismatch
  15. [15] § Results › Exploratory analyses of melatonin, sleep, and feeding history › Gut melatonin ↔ scripts/6-melatonin.ipynb, lines 78 to 101 · score 0.53 · bowel movement, spent awake, gut melatonin, feeding, age · Report a mismatch
  16. [16] § Methods › Data processing › Other variables ↔ scripts/2-EDA.ipynb, lines 187 to 206 · score 0.52 · spent awake, stool sample, feeding history, duration, clock, sleep · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,448 lines · 57 KB · no license · 5 matches

  1. import os
  2. import numpy as np
  3. import pandas as pd
  4. import matplotlib.pyplot as plt
  5. import seaborn as sns
  6. import qiime2 as q2
  7. from skbio import OrdinationResults
  8. from joypy import joyplot
  9. from scipy.optimize import curve_fit
  10. import matplotlib.lines as mlines
  11. from sklearn.model_selection import GroupShuffleSplit, GroupKFold, LeaveOneGroupOut
  12. from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
  13. from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, accuracy_score, f1_score
  16. # mapping dictionary for vizualisation labels
  17. name_mapping = {
  18. "observed_features": "Observed features",
  19. "shannon_entropy": "Shannon entropy",
  20. "pielou_evenness": "Pielou evenness",
  21. "faith_pd": "Faith phylogenetic diversity",
  22. "age_days": "Age (days)",
  23. "sex[T.Male]": "Sex (male)",
  24. "BCQ_Attunement": "Attuned care score (BCQ)",
  25. "ASQ_Composite": "Infant behavioral \n development (ASQ)",
  26. "volatility_braycurtis": "Temporal volatility: \n Bray-Curtis dissimilarity",
  27. "volatility_jaccard": "Temporal volatility: \n Jaccard similarity index",
  28. "volatility_unweighted_unifrac": "Temporal volatility: \n Unweighted UniFrac distance",
  29. "volatility_weighted_unifrac": "Temporal volatility: \n Weighted UniFrac distance",
  30. "time_since_last_bowel_movement_in_h": "Time since last \n bowel movement (h)",
  31. "prior_time_spent_awake_in_h": "Prior time spent \n awake (h)",
  32. "time_since_last_feeding_in_h": "Time since last \n feeding (h)",
  33. "duration_of_last_sleep_in_h": "Duration of last \n sleep (h)",
  34. "R2_abundance": "Alpha diversity rhythmicity: \n cosine fit score of \n observed features",
  35. "R2_abundance_evenness": "Alpha diversity rhythmicity: \n cosine fit score of \n Shannon entropy",
  36. "R2_evenness": "Alpha diversity rhythmicity: \n cosine fit score of \n Pielou evenness",
  37. "R2_biodiversity": "Alpha diversity rhythmicity: \n cosine fit score of \n Faith's phylogenetic diversity",
  38. "std_times_between_feedings_in_h": "Std. dev. of time between \n feedings (h)",
  39. "R2_bifido": "$\it{Bifidobacterium}$ rhythmicity: \n cosine fit score of the $\it{Bifidobacterium}$ genus",
  40. "R2_veillo": "$\it{Veillonella}$ rhythmicity: \n cosine fit score of the $\it{Veillonella}$ genus",
  41. "R2_esche": "$\it{Escherichia-Shigella}$ rhythmicity: \n cosine fit score of the $\it{Escherichia-Shigella}$ genus",
  42. "R2_bacteroides": "$\it{Bacteroides}$ rhythmicity: \n cosine fit score of the $\it{Bacteroides}$ genus",
  43. "R2_clostridium": "$\it{Clostridium}$ rhythmicity: \n cosine fit score of the $\it{Clostridium}$ genus"}
  45. # labels in plots with bold formatting
  46. bold_labels = ["Observed features", "Shannon entropy", "Pielou evenness", "Faith phylogenetic diversity",
  47. "Temporal volatility: \n Bray-Curtis dissimilarity", "Temporal volatility: \n Jaccard similarity index",
  48. "Temporal volatility: \n Unweighted UniFrac distance", "Temporal volatility: \n Weighted UniFrac distance",
  49. "Alpha diversity rhythmicity: \n cosine fit score of \n observed features",
  50. "Alpha diversity rhythmicity: \n cosine fit score of \n Shannon entropy",
  51. "Alpha diversity rhythmicity: \n cosine fit score of \n Pielou evenness",
  52. "Alpha diversity rhythmicity: \n cosine fit score of \n Faith's phylogenetic diversity",
  53. "$\it{Bifidobacterium}$ rhythmicity: \n cosine fit score of the $\it{Bifidobacterium}$ genus",
  54. "$\it{Veillonella}$ rhythmicity: \n cosine fit score of the $\it{Veillonella}$ genus",
  55. "$\it{Escherichia-Shigella}$ rhythmicity: \n cosine fit score of the $\it{Escherichia-Shigella}$ genus",
  56. "$\it{Bacteroides}$ rhythmicity: \n cosine fit score of the $\it{Bacteroides}$ genus",
  57. "$\it{Clostridium}$ rhythmicity: \n cosine fit score of the $\it{Clostridium}$ genus"]
  60. def process_qza_file(qza_file_path, stat_function):
  61. """
  62. Processes a .qza file to load and transform data into a processed dataframe.
  64. Parameters:
  65. qza_file_path (str): Path to the input .qza file.
  66. stat_function (function): A function that takes a dataframe as input
  67. and returns a pocessed dataframe.
  69. Returns:
  70. pd.DataFrame: The processed dataframe.
  71. """
  72. # Load the Qiime 2 Artifact as a dataframe
  73. final_feature_table_df = load_qza_artifact(qza_file_path)
  75. # Perform function
  76. final_feature_table_percent = stat_function(final_feature_table_df)
  78. return final_feature_table_percent
  81. def compute_stats_for_percentage_columns(df):
  82. """
  83. Compute the mean, median, and standard deviation of all columns
  84. with names including the string "percentage".
  86. :param df: Input pandas DataFrame
  87. :return: DataFrame with computed statistics (mean, median, std)
  88. """
  89. # Filter columns that include the string "percentage" in their names
  90. percentage_columns = [col for col in df.columns if "percentage" in col]
  92. # Select only the percentage columns
  93. df_percentage = df[percentage_columns]
  95. # Compute statistics
  96. stats = {
  97. "mean": df_percentage.mean(),
  98. "median": df_percentage.median(),
  99. "std": df_percentage.std()
  100. }
  102. # Convert the dictionary of statistics into a DataFrame
  103. stats_df = pd.DataFrame(stats)
  105. return stats_df
  108. def load_qza_artifact(path, case='metadata'):
  109. """
  110. Loads a .qza artifact and converts it into a pandas DataFrame.
  112. Parameters:
  113. path (str): Path to the .qza file.
  114. case (str): Type of artifact to load ('metadata' or 'ordination').
  116. Returns:
  117. pd.DataFrame: Extracted data as a DataFrame.
  118. """
  119. artifact = q2.Artifact.load(path)
  121. if case == 'metadata':
  122. df = artifact.view(q2.Metadata).to_dataframe()
  123. elif case == 'ordination':
  124. df = artifact.view(OrdinationResults)
  125. else:
  126. raise ValueError(f"Unsupported case: {case}. Choose 'metadata' or 'ordination'.")
  128. return df
  131. def process_pcoa(pcs, metadata):
  132. """
  133. Extracts the first three PCoA axes and merges with metadata.
  135. Parameters:
  136. pcs (pd.DataFrame): The PCoA results containing samples.
  137. metadata (pd.DataFrame): The metadata to merge.
  139. Returns:
  140. pd.DataFrame: Processed DataFrame with the first 3 PCoA axes and metadata.
  141. """
  142. pcs_processed = pcs.samples.iloc[:, :3] # Take first 3 columns
  143. pcs_processed.columns = ['Axis 1', 'Axis 2', 'Axis 3'] # Rename columns
  144. return pd.concat([pcs_processed, metadata], join="inner", axis=1) # Merge with metadata
  147. def load_tsv(file_path, index_col=0):
  148. """
  149. Loads a .tsv file into a Pandas DataFrame.
  151. Parameters:
  152. file_path (str): Path to the .tsv file.
  153. index_col (int or None): Column to use as the row labels. Default is 0 (first column).
  155. Returns:
  156. pd.DataFrame: The loaded DataFrame.
  157. """
  158. return pd.read_csv(file_path, sep='\t', index_col=index_col)
  161. def save_df_as_tsv(df, file_path, index_name="id"):
  162. """
  163. Renames the index, saves a DataFrame as a .tsv file, and ensures proper formatting.
  165. Parameters:
  166. df (pd.DataFrame): The DataFrame to be saved.
  167. file_path (str): The output file path, including the filename (e.g., 'output/data.tsv').
  168. index_name (str): The name to assign to the index column (default: 'id').
  170. Returns:
  171. None
  172. """
  173. df = df.rename_axis(index_name) # Set index name
  174. df.to_csv(file_path, sep='\t', index=True) # Save as .tsv
  177. def get_taxonomy(col_name, level='p'):
  178. """
  179. Extracts a specific taxonomic level from a taxonomy string.
  181. Parameters:
  182. col_name (str): The taxonomy string (e.g., 'k__Bacteria; p__Firmicutes').
  183. level (str): The prefix to search for (e.g., 'p' for phylum, 'g' for genus).
  185. Returns:
  186. str: The identified taxonomic rank or a fallback 'Unassigned' string.
  187. """
  188. # Map prefixes to names for the return message
  189. ranks = {'k': 'kingdom', 'p': 'phylum', 'c': 'class', 'o': 'order', 'f': 'family', 'g': 'genus'}
  190. target = f"{level}__"
  192. # Split the taxonomy string by semicolon
  193. parts = col_name.split(';')
  195. # Find the part that starts with the target
  196. for part in parts:
  197. clean_part = part.strip()
  198. if clean_part.startswith(target):
  199. return clean_part
  201. return f"Unassigned {ranks.get(level, 'rank')}"
  204. def get_tax_per_timepoint(tax_ft, metdata_ft, timepoint_col):
  205. """
  206. Helper function to get the relative abundance of a taxonomic level per timepoint.
  207. """
  209. # Merge the taxonomic feature table with metadata
  210. df_with_metadata = tax_ft.merge(metdata_ft[timepoint_col], left_index=True, right_index=True)
  211. # Group by timepoint and sum
  212. df_with_metadata_timepoint = df_with_metadata.groupby(timepoint_col).sum()
  213. # Convert to relative abundance
  214. df_with_metadata_timepoint = df_with_metadata_timepoint.div(df_with_metadata_timepoint.sum(axis=1), axis=0)
  215. # Check that rows sum to 1
  216. row_sums = df_with_metadata_timepoint.sum(axis=1)
  217. if not np.all(np.isclose(row_sums, 1.0)):
  218. print("Warning: Some total rel. abund. do not sum to 1")
  220. return df_with_metadata_timepoint
  223. def collapse_taxa(df, min_rel_abundance=0, top_n=10, keep_col='Unassigned'):
  224. """
  225. Collapses a taxonomy DataFrame by retaining top taxa and grouping the rest into 'Others'.
  227. Parameters:
  228. df (pd.DataFrame): The DataFrame to be collapsed (samples as rows, taxa as columns).
  229. min_rel_abundance (float): The minimum abundance threshold for a taxon to be kept.
  230. top_n (int): The number of most abundant taxa to keep if min_rel_abundance is 0.
  231. keep_col (str): The name of a specific column to preserve from the ranking logic (default: 'Unassigned').
  233. Returns:
  234. pd.DataFrame: A DataFrame containing the top taxa, the preserved column, and a new 'Others' column.
  235. """
  236. # 1. Identify columns to rank (exclude 'Unassigned' from the ranking logic)
  237. cols_to_rank = [c for c in df.columns if c != keep_col]
  239. if min_rel_abundance > 0:
  240. # Filter columns based on the minimum relative abundance threshold
  241. print("Filtering taxa based on the minimum relative abundance threshold: ", min_rel_abundance)
  242. top_taxa = df.columns[df.max() > min_rel_abundance]
  243. else:
  244. # 2. Get the top N taxa based on the sum across all rows
  245. print("Selecting top N taxa based on total abundance: ", top_n)
  246. top_taxa = df[cols_to_rank].sum().sort_values(ascending=False).head(top_n).index.tolist()
  248. # 3. Identify columns that are not in the top N and not the protected 'Unassigned' column
  249. others_cols = [c for c in df.columns if c not in top_taxa and c != keep_col]
  251. # 4. Build the new DataFrame
  252. # Start with the top taxa
  253. top_taxa_sorted = df[top_taxa].sum().sort_values(ascending=False).index.tolist()
  254. df_collapsed = df[top_taxa_sorted].copy()
  256. # Add 'Unassigned' if it exists in the original data
  257. if keep_col in df.columns:
  258. df_collapsed[keep_col] = df[keep_col]
  260. # Sum the remaining columns into 'Others'
  261. df_collapsed['Others'] = df[others_cols].sum(axis=1)
  263. df_final = df_collapsed.copy()
  265. return df_final
  268. def plot_taxa_stacked_bar(df, level_name, out_path, threshold=0, top_n=10, figsize=(5.5, 4)):
  269. """
  270. Plots a stacked bar chart of relative abundance.
  272. Parameters:
  273. - df: The prepared DataFrame (rows=timepoints, cols=taxa)
  274. - level_name: String (e.g., 'Phylum' or 'Genus') for labels/titles
  275. - out_path: Directory path to save the PDF
  276. - threshold: Minimum max-abundance to include a taxon (default 1%)
  277. - top_n: Number of top taxa to include
  278. - figsize: Figure size tuple
  279. """
  281. if level_name == "phylum":
  282. keep_col = 'Unassigned phylum'
  283. elif level_name == "genus":
  284. keep_col = 'Unassigned genus'
  285. else:
  286. keep_col = 'Unassigned'
  288. # Filter for abundance threshold
  289. plot_df = collapse_taxa(df, min_rel_abundance=threshold, top_n=top_n, keep_col=keep_col)
  291. fig, ax = plt.subplots(1, 1, figsize=figsize)
  293. # plotting
  294. ax.stackplot(
  295. plot_df.index,
  296. plot_df.T,
  297. labels=plot_df.columns,
  298. colors=plt.cm.tab20.colors)
  300. # Legend handling: reverse order to match the visual stack (top-to-bottom)
  301. handles, labels = ax.get_legend_handles_labels()
  302. labels = [label.split('__')[-1] if '__' in label else label for label in labels]
  304. ax.legend(handles[::-1], labels[::-1], title=level_name,
  305. loc='upper left', bbox_to_anchor=(1, 1))
  307. # Aesthetics
  308. ax.tick_params(which="both", direction="out")
  309. ax.set_ylabel("Relative abundance")
  310. ax.set_title(f"Relative abundance of gut bacterial {level_name} over time")
  312. fig.tight_layout()
  314. # Save and show
  315. save_filename = f"{out_path}/{level_name.lower()}_relative_abundance_over_time.pdf"
  316. fig.savefig(save_filename, dpi=300, bbox_inches='tight')
  317. plt.show()
  320. def compute_distance_to_centroid(df, new_col_name, x_col='Axis 1', y_col='Axis 2', z_col='Axis 3',
  321. centroid_x_col='Centroid x', centroid_y_col='Centroid y', centroid_z_col='Centroid z'):
  322. """
  323. Computes the Euclidean distance from each point (x, y, z) to its centroid (centroid_x, centroid_y, centroid_z).
  325. Parameters:
  326. df (pd.DataFrame): DataFrame containing point and centroid coordinates.
  327. x_col, y_col, z_col (str): Column names for the point coordinates.
  328. centroid_x_col, centroid_y_col, centroid_z_col (str): Column names for the centroid coordinates.
  330. Returns:
  331. pd.DataFrame: Original DataFrame with an additional column 'distance_to_centroid'.
  332. """
  333. df[new_col_name] = np.sqrt(
  334. (df[x_col] - df[centroid_x_col])**2 +
  335. (df[y_col] - df[centroid_y_col])**2 +
  336. (df[z_col] - df[centroid_z_col])**2
  337. )
  338. return df[new_col_name]
  341. def min_max_normalize(series, min_val, max_val, reverse=False):
  342. """Normalize a value using min-max scaling to range [0,1]. Optionally reverse the scale."""
  343. normalized = (series - min_val) / (max_val - min_val)
  344. return 1 - normalized if reverse else normalized
  347. def compute_sleep_score(dataframe, weights=None):
  348. """
  349. Compute a composite sleep score for each row in a dataframe.
  351. Parameters:
  352. - df: DataFrame with columns ['nighttime_sleep_duration_in_h', 'sleep_latency_in_h', 'bedtime_in_h', 'number_of_awakenings']
  353. - weights: Dictionary with weights for each variable (default: equal weights)
  355. Returns:
  356. - Additional column 'sleep_score'
  357. """
  358. df = dataframe.copy()
  360. sleep_dur = min_max_normalize(df['nighttime_sleep_duration_in_h'],
  361. df['nighttime_sleep_duration_in_h'].min(),
  362. df['nighttime_sleep_duration_in_h'].max())
  363. sleep_lat = min_max_normalize(df['sleep_latency_in_h'],
  364. df['sleep_latency_in_h'].min(),
  365. df['sleep_latency_in_h'].max(), reverse=True)
  366. bedtime = min_max_normalize(df['bedtime_in_h'],
  367. df['bedtime_in_h'].min(),
  368. df['bedtime_in_h'].max(), reverse=True)
  369. awakenings = min_max_normalize(df['number_of_awakenings'],
  370. df['number_of_awakenings'].min(),
  371. df['number_of_awakenings'].max(), reverse=True)
  373. if weights is None:
  374. weights = { 'nighttime_sleep_duration_in_h': 0.25,
  375. 'sleep_latency_in_h': 0.25,
  376. 'bedtime_in_h': 0.25,
  377. 'number_of_awakenings': 0.25 }
  379. # Compute weighted sum
  380. df['composite'] = (weights['nighttime_sleep_duration_in_h'] * sleep_dur +
  381. weights['sleep_latency_in_h'] * sleep_lat +
  382. weights['bedtime_in_h'] * bedtime +
  383. weights['number_of_awakenings'] * awakenings)
  385. return df['composite']
  388. def cosine_function(x, A, phi, C):
  389. """
  390. Cosine function to model rhythmicity. Here lambda/period = 24h.
  391. """
  392. return A * np.cos((2 * np.pi / 24) * x + phi) + C
  395. def fit_cosine_and_plot(data, timepoint_values, timepoints_colors, feature_column, ax, show_legend=True):
  396. """
  397. Fits a cosine function to the given feature column for each timepoint and subject.
  399. Parameters:
  400. - data: DataFrame containing the dataset.
  401. - timepoint_values: List of timepoints to analyze.
  402. - timepoints_colors: List of colors corresponding to timepoints.
  403. - feature_column: Column name to analyze (e.g., "observed_features", "shannon_entropy").
  404. - ax: Matplotlib Axes object to plot on.
  405. - show_legend: Boolean, whether to display the legend.
  407. Returns:
  408. - fit_scores_df: DataFrame containing R², RMSE, and MAE for each subject.
  409. - summary_stats: Summary statistics (mean, median) per timepoint.
  410. """
  412. # Store results
  413. subject_fit_scores = []
  414. timepoint_counts = {} # Store (infants, samples) per timepoint
  415. legend_handles = {}
  417. for i, timepoint in enumerate(timepoint_values):
  418. subset_data_time = data[data['timepoint'] == timepoint]
  419. color = timepoints_colors[i]
  421. n_infants = 0
  422. n_samples = 0
  424. for subject in subset_data_time['infant_id'].unique():
  425. subset_data = subset_data_time[subset_data_time['infant_id'] == subject]
  427. if subset_data.empty or len(subset_data) < 4:
  428. continue # Skip subjects with too few points
  430. n_infants += 1
  431. n_samples += len(subset_data)
  433. X = np.array(subset_data['hour_stool_sample'])
  434. y = np.array(subset_data[feature_column])
  436. # Initial parameter guesses
  437. # Amplitude (A): Half the range of y_data.
  438. # Phase (φ): Start at 0.
  439. # Offset (C): Mean of y_data.
  440. A_init = (np.max(y) - np.min(y)) / 2
  441. phi_init = 0
  442. C_init = np.mean(y)
  443. initial_guess = [A_init, phi_init, C_init]
  445. params, _ = curve_fit(cosine_function, X, y, p0=initial_guess, maxfev=5000,
  446. bounds=([0, -np.pi, np.min(y)], [np.inf, np.pi, np.max(y)]))
  447. A_fit, phi_fit, C_fit = params
  449. # Generate fit curve
  450. X_fit = np.linspace(0, 24, 100)
  451. y_fit = cosine_function(X_fit, A_fit, phi_fit, C_fit)
  453. scatter_handle = ax.scatter(X, y, color=color, alpha=0.6)
  454. ax.plot(X_fit, y_fit, color=color, linewidth=2)
  456. if timepoint not in legend_handles:
  457. legend_handles[timepoint] = scatter_handle
  459. # Compute fit quality metrics
  460. r2 = r2_score(y, cosine_function(X, A_fit, phi_fit, C_fit))
  461. rmse = np.sqrt(mean_squared_error(y, cosine_function(X, A_fit, phi_fit, C_fit)))
  462. mae = mean_absolute_error(y, cosine_function(X, A_fit, phi_fit, C_fit))
  464. # Store results
  465. subject_fit_scores.append({"infant_id": subject,
  466. "timepoint": timepoint,
  467. "R²": r2,
  468. "RMSE": rmse,
  469. "MAE": mae,
  470. "amplitude": A_fit,
  471. "constant": C_fit})
  473. timepoint_counts[timepoint] = (n_infants, n_samples)
  475. # Convert results to DataFrame
  476. fit_scores_df = pd.DataFrame(subject_fit_scores)
  478. # Compute summary stats per timepoint
  479. summary_stats = fit_scores_df.groupby("timepoint").agg({
  480. "R²": ["mean", "median"],
  481. "RMSE": ["mean", "median"],
  482. "MAE": ["mean", "median"],
  483. "amplitude": "median",
  484. "constant": "median"
  485. }).reset_index()
  487. legend_labels = [f"{tp} ({timepoint_counts[tp][0]} infants, {timepoint_counts[tp][1]} samples)"
  488. for tp in legend_handles.keys()]
  490. ax.set_title(f"Cosine fit of {feature_column.replace('_', ' ')}", fontsize=16)
  491. ax.set_xlabel("hour_stool_sample", fontsize=14)
  492. ax.set_ylabel(feature_column.replace("_", " ").title(), fontsize=14)
  493. ax.grid(True)
  494. ax.set_xticks(range(0, 25, 4))
  496. if show_legend:
  497. ax.legend(legend_handles.values(), legend_labels, title="", bbox_to_anchor=(1, 1), loc='upper left', fontsize=14)
  499. return fit_scores_df, summary_stats
  502. def fit_cosine_and_plot_scaled(data, timepoint_values, timepoints_colors, feature_column,
  503. ax, fit_scores, show_legend=True, median_amplitudes=None, median_constants=None):
  504. """
  505. Plots a scaled cosine fit using precomputed median amplitude and constant per timepoint.
  507. Parameters:
  508. - data: DataFrame containing the dataset.
  509. - timepoint_values: List of unique timepoints.
  510. - timepoints_colors: List of colors corresponding to timepoints.
  511. - feature_column: Column to analyze (e.g., "observed_features", "shannon_entropy").
  512. - ax: Matplotlib Axes object to plot on.
  513. - show_legend: Boolean, whether to display the legend.
  514. - median_amplitudes: Precomputed dictionary of median amplitude per timepoint.
  515. - median_constants: Precomputed dictionary of median constant per timepoint.
  517. Returns:
  518. - None
  519. """
  521. if median_amplitudes is None or median_constants is None:
  522. raise ValueError("median_amplitudes and median_constants must be provided.")
  524. timepoint_counts = {} # Store (infants, samples) per timepoint
  525. legend_handles = {}
  527. # Plot scaled version using precomputed median values
  528. for i, timepoint in enumerate(timepoint_values):
  529. subset_data_time = data[data['timepoint'] == timepoint]
  530. subset_R2_time = fit_scores[fit_scores['timepoint'] == timepoint]
  531. color = timepoints_colors[i]
  532. n_infants = 0
  533. n_samples = 0
  535. for subject in subset_data_time['infant_id'].unique():
  536. subset_data = subset_data_time[subset_data_time['infant_id'] == subject]
  537. subset_R2 = subset_R2_time[subset_R2_time['infant_id'] == subject]
  539. if subset_data.empty or len(subset_data) < 4:
  540. continue
  542. if subset_R2["R²"].iloc[0] < 0.5:
  543. continue
  545. n_infants += 1
  546. n_samples += len(subset_data)
  548. X = np.array(subset_data['hour_stool_sample'])
  549. y = np.array(subset_data[feature_column])
  551. # Initial parameter guesses:
  552. # Amplitude (A): Half the range of y_data.
  553. # Phase (φ): Start at 0.
  554. # Offset (C): Mean of y_data.
  555. A_init = (np.max(y) - np.min(y)) / 2
  556. phi_init = 0
  557. C_init = np.mean(y)
  558. initial_guess = [A_init, phi_init, C_init]
  560. params, _ = curve_fit(cosine_function, X, y, p0=initial_guess, maxfev=5000,
  561. bounds=([0, -np.pi, np.min(y)], [np.inf, np.pi, np.max(y)]))
  562. A_fit, phi_fit, C_fit = params
  564. X_fit = np.linspace(0, 24, 100)
  565. A_fit = list(median_amplitudes)[i][timepoint]
  566. C_fit = list(median_constants)[i][timepoint]
  567. y_fit = cosine_function(X_fit, A_fit, phi_fit, C_fit)
  569. scatter_handle = ax.scatter(X, y, color=color, alpha=0.6)
  570. ax.plot(X_fit, y_fit, color=color, linewidth=2)
  572. if timepoint not in legend_handles:
  573. legend_handles[timepoint] = scatter_handle
  575. timepoint_counts[timepoint] = (n_infants, n_samples)
  577. legend_labels = [f"{tp} ({timepoint_counts[tp][0]} infants, {timepoint_counts[tp][1]} samples)" for tp in legend_handles.keys()]
  579. ax.set_title(f"Cosine fit (R² > 0.5) of {feature_column.replace('_', ' ')} \n - Adjusted by Median Amplitude & Constant", fontsize=16)
  580. ax.set_xlabel("hour_stool_sample", fontsize=14)
  581. ax.set_ylabel(feature_column.replace("_", " ").title(), fontsize=14)
  582. ax.grid(True)
  583. ax.set_xticks(range(0, 25, 4))
  585. if show_legend:
  586. ax.legend(legend_handles.values(), legend_labels, title="", loc='lower center', fontsize=14)
  589. def blend_with_white(color, blend_factor=0.3):
  590. """
  591. Blends a given RGB color with white.
  592. `blend_factor=0.3` means 30% white, 70% original color.
  593. """
  594. white = np.array([1, 1, 1, 1]) # RGBA white
  595. return tuple((1 - blend_factor) * np.array(color) + blend_factor * white)
  598. def samples_scatter_plot(data_analysis, data_all, timepoints, output_file):
  599. """
  600. Creates a scatter plot of sample collection times with vertical bars for additional samples.
  602. Parameters:
  603. data_analysis (DataFrame): Data containing 'hour_stool_sample', 'timepoint', and 'infant_id'.
  604. data_all (DataFrame): Data containing additional samples with 'hour_stool_sample', 'timepoint', and 'infant_id'.
  605. timepoints (list): List of timepoints to include in the plot (e.g., ['2 months', '4 months', '6 months']).
  606. output_file (str): Name of the file to save the plot.
  607. """
  608. num_plots = len(timepoints)
  610. fig, axes = plt.subplots(1, num_plots, figsize=(6 * num_plots, 6), sharey=True)
  612. # Sort infants by infant_id
  613. infants = sorted(data_analysis['infant_id'].unique()) # Sorting the infant IDs
  614. infants_num = list(range(1, len(data_analysis['infant_id'].unique()) + 1))
  615. y_positions = range(len(infants))
  617. for i, timepoint in enumerate(timepoints):
  618. ax = axes[i]
  619. timepoint_data_analysis = data_analysis[data_analysis['timepoint'] == timepoint]
  621. # Count valid samples for the current timepoint
  622. num_samples = timepoint_data_analysis.dropna(subset=['hour_stool_sample']).shape[0]
  624. # Add vertical bars for the samples from `data_all` first, to be behind the scatter points
  625. timepoint_data_all = data_all[data_all['timepoint'] == timepoint]
  626. for j, infant in enumerate(infants):
  627. infant_data_all = timepoint_data_all[timepoint_data_all['infant_id'] == infant]
  628. valid_data_all = infant_data_all.dropna(subset=['hour_stool_sample'])
  630. # Increase the length of the vertical bars (set y-axis range from j-0.5 to j+0.5)
  631. ax.vlines(valid_data_all['hour_stool_sample'], j - 0.5, j + 0.5, color='black', alpha=0.4)
  633. # Now plot the scatter plot for the data_analysis (dots will be in front)
  634. for j, infant in enumerate(infants):
  635. infant_data_analysis = timepoint_data_analysis[timepoint_data_analysis['infant_id'] == infant]
  636. valid_data_analysis = infant_data_analysis.dropna(subset=['hour_stool_sample'])
  637. ax.scatter(valid_data_analysis['hour_stool_sample'], [j] * len(valid_data_analysis), alpha=1)
  639. ax.set_xticks(range(0, 24, 4)) # Add x-axis ticks for time of day
  640. ax.set_xlim(0, 24) # Set x-axis limit from 0 to 24
  641. ax.set_xlabel('Time of day (hours)', fontsize=14)
  642. ax.set_title(f'Sample collection times at {timepoint} ({num_samples} samples)', fontsize=14)
  643. ax.grid(axis='y')
  645. # Set common y-axis label
  646. axes[0].set_yticks(y_positions)
  647. axes[0].set_yticklabels(infants_num)
  648. axes[0].set_ylabel('Infants', fontsize=14)
  650. plt.tight_layout()
  652. # Save and show the plot
  653. plt.savefig(output_file, dpi=300, bbox_inches='tight')
  654. plt.show()
  657. def density_plot(data_analysis, timepoints, colors, output_file, alpha=1, ticklabel=14, labelsize=16, titlesize=16):
  658. """
  659. Creates a density plot (joyplot) for sample collection times with sample counts included in the labels.
  661. Parameters:
  662. data_analysis (DataFrame): Data containing 'hour_stool_sample' and 'timepoint' columns.
  663. timepoints (list): List of timepoints to include in the plot (e.g., ['2 months', '4 months', '6 months']).
  664. colors (list): List of colors for the density plots, one per timepoint.
  665. alpha (float): Transparency level for the density plots.
  666. ticklabel (int): Font size for tick labels.
  667. labelsize (int): Font size for axis labels.
  668. titlesize (int): Font size for the plot title.
  669. output_file (str): Name of the file to save the plot.
  670. """
  671. # Creating a dictionary of data for each timepoint and calculating sample sizes
  672. timepoint_data = {}
  673. updated_timepoints = []
  674. for timepoint in timepoints:
  675. subset = data_analysis[data_analysis['timepoint'] == timepoint]['hour_stool_sample'].dropna()
  676. timepoint_data[timepoint] = subset
  677. n_samples = len(subset)
  678. updated_timepoints.append(f"(n={n_samples} samples)")
  680. ls_colors = colors # List of colors for each timepoint
  682. # Generate joyplots
  683. fig, ax = joyplot(
  684. data=list(timepoint_data.values()), # The list of sample hour data for each timepoint
  685. labels=updated_timepoints, # Updated labels with sample sizes
  686. overlap=0, # Control overlap between plots
  687. alpha=alpha * 0.9,
  688. ylabelsize=ticklabel,
  689. xlabelsize=ticklabel,
  690. color=ls_colors, # Color for each plot
  691. linewidth=0.8, # Line width for the outlines
  692. figsize=(6, 7), # Size of the plot
  693. )
  695. # Set x-axis limits and ticks
  696. for axis in ax:
  697. axis.set_xlim(0, 24)
  698. axis.set_xticks([0, 6, 12, 18, 24])
  700. # Set y-axis label
  701. # Customize plot appearance
  702. plt.rc("font", size=labelsize)
  703. plt.title("", fontsize=titlesize, y=0.9)
  704. plt.xlabel("Time of the day (hours)", fontsize=labelsize)
  705. plt.grid(True, axis='x')
  706. plt.savefig(output_file, dpi=300, bbox_inches='tight')
  707. plt.show()
  710. def generate_color_palette(list):
  711. """
  712. Generates a custom color palette composed of lightened shades derived from
  713. selected colors within the Spectral colormap.
  715. Parameters:
  716. list (list): A list of three integers specifying the number of shades to generate
  717. from each of the three selected base colors.
  718. """
  720. # Get colors from the Spectral colormap
  721. colors = [plt.cm.Spectral(i / float(6)) for i in range(7)]
  722. timepoints_colors = [colors[0], colors[5], colors[6]]
  724. # Function to generate shades of a color
  725. def generate_shades(base_color, num_shades, lighten=True):
  726. base_rgb = np.array(base_color[:3]) # Ignore alpha channel
  727. white_or_black = np.array([1, 1, 1]) if lighten else np.array([0, 0, 0])
  729. shades = [tuple((1 - t) * base_rgb + t * white_or_black) for t in np.linspace(0, 1, num_shades)]
  730. return shades
  732. # Generate shades
  733. shades_0 = generate_shades(timepoints_colors[0], list[0], lighten=True) # Lightening colors[0]
  734. shades_5 = generate_shades(timepoints_colors[1], list[1], lighten=True) # Lightening colors[5]
  735. shades_6 = generate_shades(timepoints_colors[2], list[2], lighten=True) # Lightening colors[6]
  737. # Combine all shades
  738. final_palette = shades_0 + shades_5 + shades_6
  740. return final_palette
  743. def plot_pcoa_infants(data, column, title_suffix, ax, var_explained, final_palette="husl"):
  744. """
  745. Function to plot a PCoA scatterplot grouped by infant_id with no specific hue order
  746. and without displaying the legend.
  748. Parameters:
  749. - data: DataFrame containing the PCoA data with 'Axis 1', 'Axis 2', and the grouping column.
  750. - column: Column to use for hue (e.g., 'infant_id').
  751. - title_suffix: String to be added to the plot title.
  752. - ax: Axis object on which to plot the graph.
  753. - var_explained: List containing the percentage of variance explained by Axis 1 and Axis 2.
  754. - final_palette: Color palette to use for the plot (default is "husl").
  755. """
  757. # Get unique colors from the scatterplot
  759. data = data.sort_values(column)
  761. unique_infants = data[column].unique()
  763. if final_palette=="husl":
  764. palette = sns.color_palette("husl", len(unique_infants))
  765. else:
  766. palette = final_palette
  768. # Create the scatterplot on the given axis
  769. scatter = sns.scatterplot(data=data,
  770. x='Axis 1',
  771. y='Axis 2',
  772. hue=column,
  773. palette=palette, # Use a diverse color palette
  774. ax=ax,
  775. legend=False) # Disable the legend
  777. # Calculate centroids for each group (infant_id in this case)
  778. centroids = data.groupby(column)[['Axis 1', 'Axis 2']].mean().reset_index()
  780. color_map = dict(zip(unique_infants, palette))
  782. # Overlay filled, round centroids for each infant_id
  783. for _, row in centroids.iterrows():
  784. ax.scatter(row['Axis 1'],
  785. row['Axis 2'],
  786. color=color_map[row[column]],
  787. s=300,
  788. alpha=0.6, # Adjust alpha for shading
  789. edgecolor='black', # Add edge for better visibility
  790. linewidth=1.2)
  792. ax.set_xlabel(f"Axis 1 ({var_explained[0]})")
  793. ax.set_ylabel(f"Axis 2 ({var_explained[1]})")
  795. # Add title and labels
  796. ax.set_title(f"PCoA with distance matrix as {title_suffix}", fontsize=15)
  797. ax.grid(True)
  800. def violinplot_with_lines(data, x_column, y_column, subject_column, title, y_column_name, palette,
  801. ax=None, fontsize_title=14, save_path=None, text=None, loc=0.05):
  802. """
  803. Creates a violin plot with overlaid boxplots and subject-level connecting lines
  804. across timepoints (or any ordered x-axis variable).
  806. Parameters
  807. ----------
  808. data : pandas.DataFrame
  809. The dataset containing all required columns.
  810. x_column : str
  811. Column name for the categorical x-axis variable (e.g., timepoint).
  812. y_column : str
  813. Column name for the numerical y-axis variable.
  814. subject_column : str
  815. Column name identifying subjects for drawing repeated-measure lines.
  816. title : str
  817. Title of the plot.
  818. y_column_name : str
  819. Y-axis label.
  820. palette : list or dict
  821. Color palette to assign to each x-axis category for violin and boxplots.
  822. ax : matplotlib.axes.Axes, optional
  823. Axis to draw on. If None, a new figure and axis are created.
  824. fontsize_title : int, optional (default = 14)
  825. Font size for the plot title. Reduced automatically when using subplots.
  826. save_path : str, optional
  827. If provided, saves the figure to this path (only if a new figure is created).
  828. text : str, optional
  829. Additional annotation text to display in the top-left corner of the plot.
  830. loc : float, optional (default = 0.05)
  831. Horizontal position (axis fraction) of the annotation text.
  832. """
  834. # Sort data by x-axis values
  835. data = data.sort_values([x_column])
  836. timepoint_order = data[x_column].unique()
  838. # Get count of each timepoint
  839. count_infants = data[x_column].dropna().value_counts().sort_index()
  840. x_labels = [f"{tp} (n={count_infants[tp]})" for tp in count_infants.index]
  842. # If ax is not provided, create a new figure and axis
  843. created_fig = False
  844. if ax is None:
  845. fig, ax = plt.subplots(figsize=(8, 6))
  846. created_fig = True # Track if we created a new figure
  847. else:
  848. fontsize_title = 16 # Reduce title size for subplots
  850. # Create violin plot
  851. sns.violinplot(x=data[x_column], y=data[y_column], palette=palette, ax=ax, showfliers=False, cut=0)
  853. # Overlay boxplot
  854. sns.boxplot(x=data[x_column], y=data[y_column], saturation=0.3, width=0.4, palette=palette,
  855. boxprops={'zorder': 2}, ax=ax)
  857. # Overlay individual subject lines
  858. x_positions = {tp: i for i, tp in enumerate(timepoint_order)}
  859. for subject_id, sub_df in data.groupby(subject_column):
  860. sub_df = sub_df.sort_values(x_column)
  861. ax.plot([x_positions[tp] for tp in sub_df[x_column]], sub_df[y_column], linestyle='-',
  862. alpha=0.2, color='black')
  864. # Set titles and labels
  865. ax.set_title(title, fontsize=fontsize_title)
  866. ax.set_xlabel(None)
  867. ax.set_ylabel(y_column_name, fontsize=14)
  868. ax.grid(True)
  870. # Rotate x-axis labels if using subplots
  871. rotation_angle = 45 if not created_fig else 0
  872. ax.set_xticklabels(x_labels, fontsize=12, rotation=rotation_angle, ha="right" if rotation_angle else "center")
  874. # Add optional text annotation
  875. if text is not None:
  876. ax.text(loc, 0.95, text, transform=ax.transAxes, fontsize=9, verticalalignment='top',
  877. horizontalalignment='left', bbox=dict(facecolor='white', alpha=0.5))
  879. # Save the figure if path is provided
  880. if save_path and created_fig:
  881. plt.savefig(save_path, dpi=300, bbox_inches='tight')
  883. # If a new figure was created, show it
  884. if created_fig:
  885. plt.show()
  887. return ax # Return axis for further modifications if needed
  890. def violinplot_no_lines(data, x_column, y_column, title, y_column_name, palette,
  891. ax=None, fontsize_title=16, fontsize_axes=14, save_path=None, text=None, loc=0.05):
  892. """
  893. Creates a violin plot with overlaid boxplots.
  895. Parameters
  896. ----------
  897. data : pandas.DataFrame
  898. The dataset containing all required columns.
  899. x_column : str
  900. Column name for the categorical x-axis variable (e.g., timepoint).
  901. y_column : str
  902. Column name for the numerical y-axis variable.
  903. title : str
  904. Title of the plot.
  905. y_column_name : str
  906. Y-axis label.
  907. palette : list or dict
  908. Color palette to assign to each x-axis category for violin and boxplots.
  909. ax : matplotlib.axes.Axes, optional
  910. Axis to draw on. If None, a new figure and axis are created.
  911. fontsize_title : int, optional (default = 16)
  912. Font size for the plot title. Reduced automatically when using subplots.
  913. fontsize_axes : int, optional (default = 14)
  914. Font size for axis labels and tick labels.
  915. save_path : str, optional
  916. If provided, saves the figure to this path (only if a new figure is created).
  917. text : str, optional
  918. Additional annotation text to display in the top-left corner of the plot.
  919. loc : float, optional (default = 0.05)
  920. Horizontal position (axis fraction) of the annotation text.
  921. """
  923. # Sort data by x-axis values
  924. data = data.sort_values([x_column])
  926. # Get count of each timepoint
  927. count_samples = data[x_column].dropna().value_counts().sort_index()
  928. x_labels = [f"{tp} (n={count_samples[tp]})" for tp in count_samples.index]
  930. # If ax is not provided, create a new figure and axis
  931. created_fig = False
  932. if ax is None:
  933. fig, ax = plt.subplots(figsize=(8, 6))
  934. created_fig = True # Track if we created a new figure
  935. else:
  936. fontsize_title = 16 # Reduce title size for subplots
  938. # Create violin plot
  939. sns.violinplot(x=data[x_column], y=data[y_column], palette=palette, ax=ax, cut=0)
  941. # Overlay boxplot
  942. sns.boxplot(x=data[x_column], y=data[y_column], saturation=0.3, width=0.4, palette=palette,
  943. boxprops={'zorder': 2}, ax=ax)
  945. # Set titles and labels
  946. ax.set_title(title, fontsize=fontsize_title)
  947. ax.set_xlabel(None)
  948. ax.set_ylabel(y_column_name, fontsize=fontsize_axes)
  949. ax.grid(True)
  951. # Rotate x-axis labels if using subplots
  952. rotation_angle = 45 if not created_fig else 0
  953. ax.set_xticklabels(x_labels, fontsize=fontsize_axes, rotation=rotation_angle, ha="right" if rotation_angle else "center")
  955. # Add optional text annotation
  956. if text is not None:
  957. ax.text(loc, 0.95, text, transform=ax.transAxes, fontsize=10, verticalalignment='top',
  958. horizontalalignment='left', bbox=dict(facecolor='white', alpha=0.5))
  960. # Save the figure if path is provided
  961. if save_path and created_fig:
  962. plt.savefig(save_path, dpi=300, bbox_inches='tight')
  964. # If a new figure was created, show it
  965. if created_fig:
  966. plt.show()
  968. return ax # Return axis for further modifications if needed
  972. def scatterplot(data, x_column, y_column, timepoint_column, title, y_column_name, x_column_name, palette,
  973. ax=None, fontsize_title=16, fontsize_axes=14, save_path=None,
  974. loc_position = 'lower right', log=False):
  975. """
  976. Creates a scatter plot with regression lines and confidence intervals for different subjects or timepoints. Also annotates the legend with sample sizes.
  978. Parameters
  979. ----------
  980. data : pandas.DataFrame
  981. The dataset containing values for the x-axis, y-axis, and groups defined by `timepoint_column`.
  982. x_column : str
  983. Name of the column in `data` to be plotted on the x-axis.
  984. y_column : str
  985. Name of the column in `data` to be plotted on the y-axis.
  986. timepoint_column : str
  987. Name of the column that defines groups (e.g., timepoints) used for coloring and legend entries.
  988. title : str
  989. Title of the plot.
  990. y_column_name : str
  991. Display name (label) for the y-axis.
  992. x_column_name : str
  993. Display name (label) for the x-axis.
  994. palette : list or dict
  995. Color palette used for different timepoints. Should have at least as many entries as unique groups in `timepoint_column`.
  996. ax : matplotlib.axes.Axes, optional
  997. A Matplotlib Axes object to draw the plot on. If None, a new figure and axes are created.
  998. fontsize_title : int, optional (default=16)
  999. Font size for the plot title.
  1000. fontsize_axes : int, optional (default=14)
  1001. Font size for the axis labels.
  1002. save_path : str or None, optional
  1003. If provided, the figure will be saved to this file path.
  1004. loc_position : str or None, optional (default='lower right')
  1005. Position of the legend. If None, the legend is omitted.
  1006. """
  1008. data = data.sort_values([timepoint_column])
  1010. if ax is None:
  1011. fig, ax = plt.subplots(figsize=(5, 6))
  1013. # Scatter plot
  1014. sns.scatterplot(data=data, x=x_column, y=y_column, hue=timepoint_column, palette=palette, ax=ax)
  1016. # Add regression lines with confidence intervals for each timepoint
  1017. legend = []
  1018. updated_handles = [] # Store legend handles
  1019. i=0
  1020. for timepoint, subset in data.groupby(timepoint_column):
  1021. n_samples = len(subset)
  1022. legend.append(f"{timepoint} (n={n_samples})") # Add sample size to label
  1023. updated_handles.append(mlines.Line2D([], [], color=palette[i], marker='o', linestyle='None', markersize=8))
  1024. i+=1
  1026. sns.regplot(data=data, x=x_column, y=y_column, scatter=False, color='black', ci=95, ax=ax)
  1028. # Labels and title
  1029. ax.set_ylabel(y_column_name, fontsize=fontsize_axes)
  1030. ax.set_title(title, fontsize=fontsize_title)
  1031. ax.grid(True)
  1032. ax.set_xlabel(x_column_name, fontsize=fontsize_axes)
  1034. # Adjust legend
  1035. ax.legend_.remove()
  1036. if loc_position != None:
  1037. loc_position = loc_position
  1038. ax.legend(updated_handles, legend, title='', loc=loc_position, fontsize=12, title_fontsize=13)
  1040. # Save figure if a path is provided
  1041. if save_path:
  1042. plt.savefig(save_path, bbox_inches='tight', dpi=300)
  1044. return ax
  1047. def scatterplot_age(data, x_column, y_column, timepoint_column, title, y_column_name, x_column_name, palette,
  1048. ax=None, fontsize_title=16, fontsize_axes=14, save_path=None, loc_position='lower right'):
  1049. """
  1050. Creates a scatter plot with regression lines and confidence intervals for each group in timepoint_column. Also annotates the legend with sample sizes.
  1052. Parameters
  1053. ----------
  1054. data : pandas.DataFrame
  1055. The dataset containing values for the x-axis, y-axis, and groups defined by `timepoint_column`.
  1056. x_column : str
  1057. Name of the column in `data` to be plotted on the x-axis.
  1058. y_column : str
  1059. Name of the column in `data` to be plotted on the y-axis.
  1060. timepoint_column : str
  1061. Name of the column that defines groups (e.g., timepoints) used for coloring and legend entries.
  1062. title : str
  1063. Title of the plot.
  1064. y_column_name : str
  1065. Display name (label) for the y-axis.
  1066. x_column_name : str
  1067. Display name (label) for the x-axis.
  1068. palette : list or dict
  1069. Color palette used for different timepoints. Should have at least as many entries as unique groups in `timepoint_column`.
  1070. ax : matplotlib.axes.Axes, optional
  1071. A Matplotlib Axes object to draw the plot on. If None, a new figure and axes are created.
  1072. fontsize_title : int, optional (default=16)
  1073. Font size for the plot title.
  1074. fontsize_axes : int, optional (default=14)
  1075. Font size for the axis labels.
  1076. save_path : str or None, optional
  1077. If provided, the figure will be saved to this file path.
  1078. loc_position : str or None, optional (default='lower right')
  1079. Position of the legend. If None, the legend is omitted.
  1080. """
  1082. data = data.sort_values([timepoint_column])
  1084. if ax is None:
  1085. fig, ax = plt.subplots(figsize=(8, 6))
  1087. # Scatter plot
  1088. sns.scatterplot(data=data, x=x_column, y=y_column, hue=timepoint_column, palette=palette, ax=ax)
  1090. # Regression lines per timepoint
  1091. legend = []
  1092. updated_handles = []
  1094. for i, (timepoint, subset) in enumerate(data.groupby(timepoint_column)):
  1095. n_samples = len(subset)
  1096. legend.append(f"{timepoint} (n={n_samples})")
  1097. updated_handles.append(mlines.Line2D([], [], color=palette[i], marker='o', linestyle='None', markersize=8))
  1099. sns.regplot(
  1100. data=subset,
  1101. x=x_column,
  1102. y=y_column,
  1103. scatter=False,
  1104. color=palette[i],
  1105. ci=95,
  1106. ax=ax
  1107. )
  1109. # Labels and title
  1110. ax.set_xlabel(x_column_name, fontsize=fontsize_axes)
  1111. ax.set_ylabel(y_column_name, fontsize=fontsize_axes)
  1112. ax.set_title(title, fontsize=fontsize_title)
  1113. ax.grid(True)
  1115. # Adjust legend
  1116. ax.legend_.remove()
  1117. if loc_position is not None:
  1118. ax.legend(updated_handles, legend, title='', loc=loc_position, fontsize=12, title_fontsize=13)
  1120. # Save figure if a path is provided
  1121. if save_path:
  1122. plt.savefig(save_path, bbox_inches='tight', dpi=300)
  1124. return ax
  1127. def plot_estimates(models, titles, x_pos, save_path):
  1128. """
  1129. Creates a multi-panel forest plot for rhythmicity models.
  1131. Parameters:
  1132. -----------
  1133. models : list
  1134. List of fitted model objects.
  1135. titles : list of str
  1136. Titles for each subplot.
  1137. x_pos : list of float
  1138. X-axis positions for the p-value annotations in each subplot.
  1139. figures_path : str
  1140. Directory path where the PDF will be saved.
  1141. """
  1143. predictors = models[0].params.drop(['Intercept', 'Group Var'], errors='ignore').rename(index=name_mapping).index
  1145. fig, axes = plt.subplots(1, len(models), figsize=(22, len(predictors)), sharey=True)
  1147. # Ensure axes is iterable if only one model is provided
  1148. if len(models) == 1:
  1149. axes = [axes]
  1151. for i, (model, ax, title) in enumerate(zip(models, axes, titles)):
  1152. params = model.params.drop(['Intercept', 'Group Var'], errors='ignore').rename(index=name_mapping)
  1153. conf = model.conf_int().drop(['Intercept', 'Group Var'], errors='ignore').rename(index=name_mapping)
  1154. p_values = model.pvalues.drop(['Intercept', 'Group Var'], errors='ignore').rename(index=name_mapping)
  1156. # Align predictors to a consistent order
  1157. params = params.reindex(predictors)
  1158. conf = conf.reindex(predictors)
  1159. p_values = p_values.reindex(predictors)
  1161. # Compute error bars
  1162. lower = params - conf[0]
  1163. upper = conf[1] - params
  1165. # Plot
  1166. ax.errorbar(params.values,
  1167. params.index,
  1168. xerr=[lower.values, upper.values],
  1169. fmt='o',
  1170. color='black',
  1171. ecolor='gray',
  1172. capsize=3)
  1174. ax.axvline(x=0, color='blue', linestyle='--')
  1175. ax.set_title(title, fontsize=16)
  1176. ax.set_xlabel("β Estimate", fontsize=14)
  1177. ax.set_ylabel("")
  1178. ax.set_yticklabels([])
  1179. ax.grid(axis='y')
  1180. ax.set_ylim(-0.5, len(predictors) - 0.5)
  1181. for tick in ax.get_yticklabels():
  1182. if tick.get_text() in bold_labels:
  1183. tick.set_fontweight('bold')
  1185. # Add p-values as text annotations next to the points
  1186. pos = x_pos[i]
  1187. for j, (coef, p_val) in enumerate(zip(params.values, p_values.values)):
  1188. if np.isnan(p_val):
  1189. label = "p = nan"
  1190. else:
  1191. label = f"p = {p_val:.3f}" if p_val >= 0.001 else "p < 0.001"
  1192. ax.text(pos, j + 0.1, label, verticalalignment='center', fontsize=10)
  1194. # Label only the first axis
  1195. axes[0].set_yticks(range(len(predictors)))
  1196. axes[0].set_yticklabels(predictors, fontsize=14)
  1197. axes[0].set_ylabel('Predictors', fontsize=14)
  1199. plt.savefig(f"{save_path}", dpi=300, bbox_inches='tight')
  1200. plt.show()
  1203. def plot_estimates_vertically(models, titles, save_path):
  1204. """
  1205. Creates a multi-panel forest plot for rhythmicity models.
  1207. Parameters:
  1208. -----------
  1209. models : list
  1210. List of fitted model objects.
  1211. titles : list of str
  1212. Titles for each subplot.
  1213. save_path : str
  1214. Directory path where the PDF will be saved.
  1215. """
  1217. n_models = len(models)
  1218. if n_models <= 2:
  1219. x_size = 5
  1220. else:
  1221. x_size = 12
  1223. fig, axes = plt.subplots(n_models, 1, figsize=(x_size, 5 * n_models), sharex=True)
  1225. # Ensure axes is iterable if only one model is provided
  1226. if n_models == 1:
  1227. axes = [axes]
  1228. x_pos = 50
  1229. else:
  1230. x_pos = 0.3
  1232. for i, (model, ax, title) in enumerate(zip(models, axes, titles)):
  1233. params = model.params.drop(['Intercept', 'Group Var'], errors='ignore').rename(index=name_mapping)
  1234. conf = model.conf_int().drop(['Intercept', 'Group Var'], errors='ignore').rename(index=name_mapping)
  1235. p_values = model.pvalues.drop(['Intercept', 'Group Var'], errors='ignore').rename(index=name_mapping)
  1237. # Compute error bars
  1238. lower = params - conf[0]
  1239. upper = conf[1] - params
  1241. # Plot
  1242. ax.errorbar(params.values,
  1243. params.index,
  1244. xerr=[lower.values, upper.values],
  1245. fmt='o',
  1246. color='black',
  1247. ecolor='gray',
  1248. capsize=3)
  1250. ax.axvline(x=0, color='blue', linestyle='--')
  1251. ax.set_title(title, fontsize=16, pad=10)
  1252. ax.set_xlabel("β Estimate", fontsize=14)
  1253. ax.set_ylabel("Predictors", fontsize=14)
  1254. ax.grid(axis='y')
  1256. # Show predictor names for each model
  1257. ax.set_yticks(range(len(params.index)))
  1258. ax.set_yticklabels(params.index, fontsize=12)
  1259. ax.set_ylim(-0.5, len(params.index) - 0.5)
  1261. for tick in ax.get_yticklabels():
  1262. if tick.get_text() in bold_labels:
  1263. tick.set_fontweight('bold')
  1265. # Add p-values as text annotations next to the points
  1266. for j, (coef, p_val) in enumerate(zip(params.values, p_values.values)):
  1267. ax.text(x_pos, j + 0.2, f"p = {p_val:.3f}" if p_val >= 0.001 else "p < 0.001",
  1268. verticalalignment='center', fontsize=10)
  1270. # Adjust spacing
  1271. plt.subplots_adjust(hspace=0.5)
  1273. # Save and show using the passed save_path
  1274. plt.savefig(save_path, dpi=300, bbox_inches='tight')
  1275. plt.show()
  1278. def prepare_features(feature_table, data_analysis, data_timepoint_sleep, column='babySQUID'):
  1279. """
  1280. Prepare feature and label arrays for model training by merging feature data,
  1281. metadata, and target sleep-quality measurements.
  1282. This function performs the following steps:
  1283. 1. Merges `feature_table` with `data_analysis` to add `infant_id` and `timepoint`
  1284. based on their shared index.
  1285. 2. Groups features by `infant_id` and `timepoint` and computes the mean for each group.
  1286. 3. Merges the aggregated features with sleep-quality data (`data_timepoint_sleep`)
  1287. using the specified target column.
  1288. 4. Extracts feature matrix `X`, target vector `y`, and group labels for grouped
  1289. cross-validation.
  1291. Parameters
  1292. ----------
  1293. feature_table : pandas.DataFrame
  1294. A dataframe containing extracted features with an index matching that of
  1295. `data_analysis`.
  1296. data_analysis : pandas.DataFrame
  1297. Contains at least the columns `infant_id` and `timepoint`, aligned with
  1298. `feature_table` via index.
  1299. data_timepoint_sleep : pandas.DataFrame
  1300. Table containing sleep-related outcomes, including the specified `column`
  1301. (default: 'babySQUID'), along with `infant_id` and `timepoint`.
  1302. column : str, default='babySQUID'
  1303. The name of the target column to extract from `data_timepoint_sleep`.
  1304. """
  1305. # Merge feature_table with data_analysis on index and specified columns
  1306. feature_table_temp = feature_table.merge(
  1307. data_analysis[['infant_id', 'timepoint']],
  1308. left_index=True,
  1309. right_index=True
  1310. )
  1312. # Group by infant_id and timepoint, then compute mean
  1313. feature_table_mean = feature_table_temp.groupby(['infant_id', 'timepoint']).mean()
  1315. # Merge with sleep quality data
  1316. feature_table_all = feature_table_mean.merge(
  1317. data_timepoint_sleep[['infant_id', 'timepoint', column]],
  1318. left_index=True,
  1319. right_on=['infant_id', 'timepoint'],
  1320. how='inner'
  1321. )
  1323. # Set and reset index
  1324. feature_table_all.set_index(['infant_id', 'timepoint'], inplace=True)
  1325. feature_table_all.reset_index(inplace=True)
  1327. # Extract features and labels
  1328. features = feature_table_all.drop(columns=[column, 'infant_id', 'timepoint'])
  1329. X = features.values
  1330. y = feature_table_all[column].values
  1331. groups = feature_table_all['infant_id']
  1333. return X, y, groups, features
  1336. def evaluate_model(X, y, groups, train_test_split, task_type="regression", n_splits=5, test_size=None):
  1337. """
  1338. Evaluate a machine learning model using grouped cross-validation.
  1340. This function supports both regression and classification tasks and allows for
  1341. multiple grouped splitting strategies (GroupShuffleSplit, GroupKFold,
  1342. LeaveOneGroupOut). It trains a Random Forest model on each split, computes the
  1343. appropriate evaluation metrics, and returns the best-performing model based on
  1344. R² (regression) or accuracy (classification).
  1346. Parameters
  1347. ----------
  1348. X : array-like of shape (n_samples, n_features)
  1349. Feature matrix.
  1350. y : array-like of shape (n_samples,)
  1351. Target vector.
  1352. groups : array-like of shape (n_samples,)
  1353. Group labels used for ensuring train/test splits respect group boundaries.
  1354. train_test_split : {"GroupShuffleSplit", "GroupKFold", "LeaveOneGroupOut"}
  1355. The cross-validation strategy to apply.
  1356. task_type : {"regression", "classification"}, default="regression"
  1357. Determines which model and metrics to use.
  1358. n_splits : int, default=5
  1359. Number of folds or splits (ignored for LeaveOneGroupOut).
  1360. test_size : float, optional
  1361. Test set proportion for GroupShuffleSplit. Ignored for other strategies.
  1362. """
  1363. # Select model based on task
  1364. if task_type == "regression":
  1365. model = RandomForestRegressor(n_estimators=500, random_state=42)
  1366. elif task_type == "classification":
  1367. model = RandomForestClassifier(n_estimators=500, random_state=42)
  1368. else:
  1369. raise ValueError(f"Unsupported task type: {task_type}")
  1371. # Choose splitting strategy
  1372. match train_test_split:
  1373. case 'GroupShuffleSplit':
  1374. random_splits = GroupShuffleSplit(n_splits=n_splits, test_size=test_size, random_state=42)
  1375. case 'GroupKFold':
  1376. random_splits = GroupKFold(n_splits=n_splits)
  1377. case 'LeaveOneGroupOut':
  1378. random_splits = LeaveOneGroupOut()
  1379. case _:
  1380. raise ValueError(f"Unsupported split type: {train_test_split}")
  1382. # Score containers
  1383. r2_scores = []
  1384. mse_scores = []
  1385. mae_scores = []
  1386. acc_scores = []
  1387. f1_scores = []
  1389. best_score = -np.inf
  1390. best_model = None
  1391. best_y_pred = None
  1392. best_y_test = None
  1393. best_X_test = None
  1395. for i, (train_index, test_index) in enumerate(random_splits.split(X, y, groups)):
  1396. X_train, X_test = X[train_index], X[test_index]
  1397. y_train, y_test = y[train_index], y[test_index]
  1399. model.fit(X_train, y_train)
  1400. y_pred = model.predict(X_test)
  1402. if task_type == "regression":
  1403. if len(y_test) < 2:
  1404. r2 = np.nan
  1405. else:
  1406. r2 = r2_score(y_test, y_pred)
  1407. r2_scores.append(r2)
  1409. mse = mean_squared_error(y_test, y_pred)
  1410. mae = mean_absolute_error(y_test, y_pred)
  1411. mse_scores.append(mse)
  1412. mae_scores.append(mae)
  1414. current_score = r2 if not np.isnan(r2) else -np.inf
  1416. elif task_type == "classification":
  1417. acc = accuracy_score(y_test, y_pred)
  1418. f1 = f1_score(y_test, y_pred, average='weighted')
  1419. acc_scores.append(acc)
  1420. f1_scores.append(f1)
  1421. current_score = acc # Or use f1
  1423. # Save best model
  1424. if current_score > best_score:
  1425. best_score = current_score
  1426. best_model = model
  1427. best_y_pred = y_pred
  1428. best_y_test = y_test
  1429. best_X_test = X_test
  1431. # Summary
  1432. print("\n--- Cross-Validation Summary ---")
  1433. if task_type == "regression":
  1434. avg_r2 = np.nanmean(r2_scores)
  1435. avg_mse = np.mean(mse_scores)
  1436. avg_mae = np.mean(mae_scores)
  1437. print(f"Average R²: {avg_r2:.3f}")
  1438. print(f"Average MSE: {avg_mse:.3f}")
  1439. print(f"Average MAE: {avg_mae:.3f}")
  1440. else:
  1441. avg_acc = np.mean(acc_scores)
  1442. std_acc = np.std(acc_scores)
  1443. avg_f1 = np.mean(f1_scores)
  1444. std_f1 = np.std(f1_scores)
  1445. print(f"Average Accuracy: {avg_acc:.3f} (±{std_acc:.4f})")
  1446. print(f"Average F1 Score (weighted): {avg_f1:.3f} (±{std_f1:.4f})\n")
  1448. return best_model, best_y_pred, best_y_test, best_X_test

functions_script.py at commit bbecc5f, no license · at the source

Overview

Authors: Fannie Kerff1, Christophe Mühlematter2, Anja Adamov1, Deborah Fast2, Serafina Plüss1, Petra Zimmermann3,4,5, Salome Kurth2,6, Nicholas A Bokulich1
ORCID iDs: Fannie Kerff, Christophe Mühlematter, Anja Adamov, Petra Zimmermann, Salome Kurth, Nicholas A Bokulich
  1. Department of Health Sciences and Technology, ETH Zurich, Zurich, Switzerland
  2. Department of Psychology, University of Fribourg, Fribourg, Switzerland
  3. Faculty of Health Science and Medicine, University of Lucerne, Lucerne, Switzerland
  4. Department of Paediatrics, Children's Hospital of Central Switzerland, Lucerne, Switzerland
  5. Department of Community Health, Faculty of Science and Medicine, University of Fribourg, Fribourg, Switzerland
  6. Department of Pulmonology, University Hospital Zurich, Zurich, Switzerland
Institutions: ETH Zurich (Switzerland); University of Fribourg (Switzerland); University of Lucerne (Switzerland); University Hospital Zurich (Switzerland)
Journal: Gut microbes reports, volume 3, issue 1, article 2649096
Dates: received 12 August 2025; accepted 17 March 2026; published online 29 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1080/29933935.2026.2649096 · PMID 41924305 · PMCID PMC13037436 · OpenAlex W7143401029
Open access: diamond, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: human (organism), developmental (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning
Keywords: Gut microbiota, infant sleep, microbiota–gut–brain axis, melatonin, diurnal rhythmicity
Topic: Neuroscience of respiration and sleep (Endocrine and Autonomic Systems, Neuroscience), according to OpenAlex
Funding: Swiss National Science Foundation (181279)
Citations: cited by 2 papers (Europe PMC); 73 references in the paper

Abstract

Emerging evidence supports a bidirectional relationship between the gut microbiome and sleep, which is partly mediated by the microbiota‒gut‒brain axis. Infancy is a critical window for the establishment of both the gut microbiome and sleep regulation, which we hypothesize to be linked across both short (diurnal) and long (monthly) time scales. In this longitudinal study, we investigated associations between gut microbiota development and sleep patterns in 20 infants at 2, 4, and 6 months of age (n = 163 samples). Infants were continuously monitored across 48-h sampling periods. The gut microbiota profiles were characterized using 16S rRNA gene sequencing; gut melatonin concentrations were measured; sleep data were collected via wearable actimetry, 24-h parent-reported sleep diaries, and the Brief Infant Sleep Questionnaire; and parenting style and behavioral development were assessed. In some infants, bacterial diversity followed diurnal rhythmic patterns. While bacterial rhythmicity was not significantly associated with sleep rhythmicity, infants with higher microbial alpha diversity showed more robust sleep patterns. Infant age emerged as the strongest predictor of gut microbial diversity and melatonin levels. Our findings suggest that gut microbial establishment may support the maturation of sleep‒wake rhythms in early infancy. Further research is needed to elucidate mechanistic roles of the gut microbiome in sleep development.

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

Repository

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

bokulich-publications/microbiome-sleep-rhythmicity

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: bbecc5f41408406edcd5b4d3e7864856c153c18a, 15 March 2026
Languages: Shell (19), Jupyter (9), Python (1)
Size: 43 files, 29 scripts
Software Heritage: not archived
Found in: “Data availability statement”
Holds: README, environment (environment.yml), 9 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (9 files), pandas (5 files), statsmodels (5 files), NumPy (3 files), scikit-learn (2 files), SciPy (2 files), seaborn (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
30 files
At the source: github.com/bokulich-publications/microbiome-sleep-rhythmicity

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

Tracing map

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

What the map holds:

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

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

Data

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

Data availability statement

Data and code availability

Newly generated sequencing data (16S rRNA) from stool samples in this study are publicly available in the NCBI Sequence Read Archive (SRA) under accession code PRJEB104111.

All original code supporting the findings of this study are publicly available on GitHub at: https://github.com/bokulich-publications/microbiome-sleep-rhythmicity

Any additional information required to reanalyze the data reported in this paper is available from the corresponding authors upon request.

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, 8 authors, 5 keywords, 1 funder, 70 references.

Cite

This paper

Kerff, F., Mühlematter, C., Adamov, A., Fast, D., Plüss, S., Zimmermann, P., Kurth, S., & Bokulich, N. A. (2026). The gut microbiota and sleep in infants: a focus on diurnal rhythmicity patterns. Gut microbes reports, 3(1), 2649096. https://doi.org/10.1080/29933935.2026.2649096

BibTeX

@article{kerff2026gut,
author = {Kerff, Fannie and Mühlematter, Christophe and Adamov, Anja and Fast, Deborah and Plüss, Serafina and Zimmermann, Petra and Kurth, Salome and Bokulich, Nicholas A},
title = {{The gut microbiota and sleep in infants: a focus on diurnal rhythmicity patterns}},
journal = {Gut microbes reports},
year = {2026},
month = mar,
volume = {3},
number = {1},
pages = {2649096},
publisher = {Taylor \& Francis},
issn = {2993-3935},
doi = {10.1080/29933935.2026.2649096},
url = {https://doi.org/10.1080/29933935.2026.2649096},
pmid = {41924305},
pmcid = {PMC13037436}
}

RIS

TY - JOUR
AU - Kerff, Fannie
AU - Mühlematter, Christophe
AU - Adamov, Anja
AU - Fast, Deborah
AU - Plüss, Serafina
AU - Zimmermann, Petra
AU - Kurth, Salome
AU - Bokulich, Nicholas A
TI - The gut microbiota and sleep in infants: a focus on diurnal rhythmicity patterns
T2 - Gut microbes reports
J2 - Gut Microbes Rep
PY - 2026
DA - 2026/03/29
VL - 3
IS - 1
SP - 2649096
SN - 2993-3935
PB - Taylor & Francis
DO - 10.1080/29933935.2026.2649096
UR - https://doi.org/10.1080/29933935.2026.2649096
LA - en
ER -

CSL-JSON

{
"id": "10.1080/29933935.2026.2649096",
"type": "article-journal",
"title": "The gut microbiota and sleep in infants: a focus on diurnal rhythmicity patterns",
"container-title": "Gut microbes reports",
"author": [
{
"family": "Kerff",
"given": "Fannie"
},
{
"family": "Mühlematter",
"given": "Christophe"
},
{
"family": "Adamov",
"given": "Anja"
},
{
"family": "Fast",
"given": "Deborah"
},
{
"family": "Plüss",
"given": "Serafina"
},
{
"family": "Zimmermann",
"given": "Petra"
},
{
"family": "Kurth",
"given": "Salome"
},
{
"family": "Bokulich",
"given": "Nicholas A"
}
],
"container-title-short": "Gut Microbes Rep",
"volume": "3",
"issue": "1",
"page": "2649096",
"DOI": "10.1080/29933935.2026.2649096",
"PMID": "41924305",
"PMCID": "PMC13037436",
"ISSN": "2993-3935",
"publisher": "Taylor & Francis",
"URL": "https://doi.org/10.1080/29933935.2026.2649096",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
29
]
]
}
}

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.1101/2025.07.11.664337
The gut microbiota and sleep in infants: a focus on diurnal rhythmicity patterns
Journal: bioRxiv (preprint)
In common: developmental, 51 references, 4 authors
[2] doi:10.1371/journal.pone.0345651 [code]
Non-concussive head impacts sustained during American football correlate with changes in gut microbiome diversity and composition.
Journal: PloS one
In common: statsmodels, seaborn, scikit-learn, 4 other tools, 5 references
[3] doi:10.1038/s44323-026-00071-7 [code]
Tracing infant sleep neurophysiology longitudinally from 3 to 6 months: EEG insights into brain development.
Journal: Npj biological timing and sleep
In common: developmental, 3 references, author Salome Kurth
[4] doi:10.1186/s40168-025-02095-w [code]
Distinct gut microbiome characteristics and dynamics in patients with Parkinson's disease based on the presence of premotor rapid-eye movement sleep behavior disorders.
Journal: Microbiome
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 6 references
[5] doi:10.1016/j.xcrm.2024.101845 [code]
The neonatal gut microbiota: A role in the encephalopathy of prematurity.
Journal: Cell reports. Medicine
In common: developmental, 7 references
[6] doi:10.1126/sciadv.adj0400 [code]
A digital twin of the infant microbiome to predict neurodevelopmental deficits.
Journal: Science advances
In common: statsmodels, seaborn, scikit-learn, 4 other tools, developmental, 3 references
[7] doi:10.3389/fnut.2026.1776896
Maternal docosahexaenoic acid supplementation shapes offspring gut microbiota to modulate the gut-brain axis in a sow-piglet model.
Journal: Frontiers in nutrition
In common: developmental, 6 references
[8] doi:10.1186/s13054-023-04531-5 [code]
Modifications of lung microbiota structure in traumatic brain injury ventilated patients according to time and enteral feeding formulas: a prospective randomized study.
Journal: Critical care (London, England)
In common: statsmodels, seaborn, pandas, 3 other tools, 3 references
[9] doi:10.1038/s41598-023-44676-5
Superior performance of biofilm versus planktonic Limosilactobacillus reuteri in protection of the intestines and brain in a piglet model of necrotizing enterocolitis.
Journal: Scientific reports
In common: developmental, 6 references
[10] doi:10.3389/fnins.2023.1212791
Physiological and cognitive changes after treatments of cyclophosphamide, methotrexate, and fluorouracil: implications of the gut microbiome and depressive-like behavior.
Journal: Frontiers in neuroscience
In common: 6 references

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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1080/29933935.2026.2649096.svg)](https://openscicode.org/badge/doi/10.1080/29933935.2026.2649096)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1080/29933935.2026.2649096.svg
   :target: https://openscicode.org/badge/doi/10.1080/29933935.2026.2649096
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1080/29933935.2026.2649096"><img src="https://openscicode.org/badge/doi/10.1080/29933935.2026.2649096.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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