OSCR

Faecalibacterium prausnitzii, depleted in the Parkinson's disease microbiome, improves motor deficits in α-synuclein overexpressing mice.

Code ↔ Paper

9 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 9 matches
  1. [1] § Methods › Motor and gastrointestinal function testing ↔ notebooks/Fp_motor_pooling_MLE_analysis.ipynb, lines 369–492 · score 0.87 · bootstrap resampling, bias corrected, bootstrap distribution, bootstrap sample, Confidence intervals, accelerated
  2. [2] § Methods › Fecal microbiome analysis with shotgun metagenomics ↔ Main-Analysis-SGB.Rmd, lines 271–304 · score 0.73 · addMDS, PCoA, genotype treatment, mia, beta, variance
  3. [3] § Methods › RNA sequencing ↔ notebook/shell_scripts/fastq_processing.sh, lines 42–80 · score 0.71 · FeatureCounts, Primary Assembly, M36, GENCODE, genome, RNA
  4. [4] § Results › F. prausnitzii supplementation remodels the gene expression profile of the large intestine ↔ notebook/bulk-rnaseq-analysis.qmd, lines 788–866 · score 0.69 · Slc8a3, Tnfrsf12a, gene expression, Edn3, Thy1 ASO, intestinal
  5. [5] § Methods › RNA sequencing ↔ notebook/R_scripts/quick_enrich.R, lines 14–69 · score 0.67 · log2 fold change, clusterProfiler, ontology, enrichment, Gene
  6. [6] § Results › F. prausnitzii supplementation remodels the gene expression profile of the large intestine ↔ notebook/bulk-rnaseq-analysis.qmd, lines 788–866 · score 0.66 · Cd8a, gene expression, CCL5, Tfrc, Thy1 ASO, large intestine
  7. [7] § Results › F. prausnitzii treatment aligns the Thy1-ASO microbiome to more closely resemble a WT profile ↔ Main-Analysis-SGB.Rmd, lines 595–651 · score 0.63 · beam traversal errors, MetaPhlAn, ASO controls, mediators, SGBs, Mediation
  8. [8] § Results › F. prausnitzii supplementation remodels the gene expression profile of the large intestine ↔ notebook/bulk-rnaseq-analysis.qmd, lines 681–785 · score 0.54 · log2 fold change, gene expression, Pearson, scatter, Box, large intestine
  9. [9] § Results › F. prausnitzii treatment aligns the Thy1-ASO microbiome to more closely resemble a WT profile ↔ Main-Analysis-SGB.Rmd, lines 489–534 · score 0.52 · bead expulsion, GI symptom, mediator, beam, mediation, errors

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

R Markdown · 961 lines · 35 KB · no license · 3 matches

  1. ---
  2. title: "*Faecalibacterium prausnitzii*, depleted in the Parkinson’s disease microbiome, improves motor deficits in alpha-synuclein overexpressing mice"
  3. subtitle: "Data: Metaphlan 4.1.1 relative abundances"
  4. author:
  5. - name: Anastasiya Moiseyenko
  6. - name: Giacomo Antonello
  7. - name: Aubrey Schonhoff
  8. - name: Kaelyn Long
  9. - name: Joseph Boktor
  10. - name: Blake Dirks
  11. - name: Anastasiya Oguienko
  12. - name: Alex Villoria-Winnett
  13. - name: Rustem Ismagilov
  14. - name: Rosa Krajmalnik-Brown
  15. - name: Nicola Segata
  16. - name: Levi Waldron
  17. - name: Sarkis K. Mazmanian
  18. date: "`r format(Sys.time(), '%B %d, %Y')`"
  19. output:
  20. rmdformats::robobook:
  21. self_contained: true
  22. toc_depth: 3
  23. toc_float: false
  24. use_bookdown: true
  25. highlight: haddock
  26. params:
  27. basic_outdir: "results"
  28. assayTested: "relabundance"
  29. microbiome_data_type: "SGB"
  30. beta_dist: "bray"
  31. beta_MDS: "MDS"
  32. ---
  34. ```{r setup, include=FALSE}
  35. knitr::opts_chunk$set(warning = FALSE, message = FALSE)
  37. # packages for tabular data manipulation
  38. library(tidyverse)
  40. # microbiome specific packages
  41. library(maaslin3) # this version needs a custom fix in a fork I made pending approval as of Jul 28 2025. install with remotes::install_github("g-antonello/maaslin3", force = TRUE)
  42. library(mia)
  43. library(biobakeryUtils)
  44. # some packages for data visualization (plots and tables)
  45. library(miaViz)
  46. library(ggpubr)
  47. library(kableExtra)
  48. library(ggtext) # enhanced text handling in ggplot
  50. # MAC's package to download data
  51. library(parkinsonsMetagenomicData)
  52. # library to clean linear models output
  53. library(broom)
  54. library(ggh4x) # for complex ggplot2 faceting
  55. # remotes::install_github("cmartin/ggConvexHull")
  56. library(ggConvexHull) # for convex hulls geometries
  57. # set ggplot2 plot themes
  58. theme_set(theme_light())
  59. # set seed for reproducibility
  60. set.seed(1234)
  61. # bugsigdb packages
  62. library(bugsigdbr)
  63. library(bugSigSimple)
  65. # utility functions for biobakery and other microbiome analyses
  66. library(biobakeryUtils)
  67. ###################################################################
  68. # pairwise PERMANOVA source code
  69. source("https://raw.githubusercontent.com/pmartinezarbizu/pairwiseAdonis/refs/heads/master/pairwiseAdonis/R/pairwise.adonis.R")
  70. ###################################################################
  72. ###################################################################
  73. # function to colorize text of ggplot object matching with a pattern
  74. # This function will apply ALL color rules to EACH label string
  75. colorize_labels <- function(labels, mapping) {
  76. # Start with the original labels
  77. colored_labels <- labels
  79. # Loop through each pattern-color rule in your map
  80. for (pattern in names(mapping)) {
  81. color <- mapping[[pattern]]
  83. # The text to display is the pattern without the escape slashes
  84. display_text <- gsub("\\\\", "", pattern)
  86. # Apply the coloring rule. This will modify labels that have already been partially colored.
  87. colored_labels <- gsub(
  88. pattern = pattern,
  89. replacement = paste0("<span style='color:", color, ";'>", display_text, "</span>"),
  90. x = colored_labels,
  91. fixed = TRUE
  92. )
  93. }
  95. return(colored_labels)
  96. }
  98. ###################################################################
  99. # load data for analyses and create output directory
  100. #####################################
  101. # Load data for analyses
  102. MetaPhlAn.tse <- read_TSE_from_dir("Data/MetaPhlAn.tse/")
  104. # create base output directory
  105. dir.create(file.path(params$basic_outdir, params$microbiome_data_type),
  106. showWarnings = FALSE, recursive = TRUE)
  108. # make a study-wide color palette for the variables tested
  109. overall_color_palette <- c("WT Control" = "#b2df8a", "Thy1-ASO Control" = "#a6cee3", "Thy1-ASO *F. prausnitzii*" = "#1f78b4", "Thy1-ASO F. prausnitzii" = "#1f78b4")
  112. # load Genome Taxonomy Data Base (originally from the MetaPhlAn GitHub release)
  113. GTDB.df <- read_tsv("Data/GTDB_202403.tsv")
  115. # Initialize figure indices
  116. SuppFigIndex <- 1
  117. MainFigIndex <- 1
  118. SuppTabIndex <- 1
  119. MainTabIndex <- 1
  121. # print input data
  122. MetaPhlAn.tse
  123. ```
  125. # Preliminary statistics
  127. ## Uniform reads distribution per condition
  129. ```{r, fig.height=7, fig.width=7}
  130. normality_test <- shapiro.test(MetaPhlAn.tse@colData$number_reads)
  132. hist_reads <- ggplot(colData(MetaPhlAn.tse), aes(x = number_reads/1e6, color = genotype_treatment, fill = genotype_treatment)) +
  133. scale_color_manual(values = overall_color_palette, aesthetics = c("fill", "color")) +
  134. geom_density(alpha = 0.1) +
  135. labs(
  136. x = "Million reads per sample retained",
  137. y = "Density",
  138. color = "Genotype and Treatment",
  139. fill = "Genotype and Treatment",
  140. caption = sprintf("Shapiro-Wilk: W = %.3f, P = %.3f", normality_test$statistic, normality_test$p.value)
  141. )
  143. test_reads <- pairwise.t.test(MetaPhlAn.tse@colData$number_reads, MetaPhlAn.tse@colData$genotype_treatment, pool.sd = FALSE,p.adjust.method = NULL) %>%
  144. broom::tidy()
  146. reads_per_condition <- ggarrange(hist_reads, ggtexttable(test_reads, rows = NULL), common.legend = TRUE, labels = "AUTO", nrow = 2, heights = c(1, 0.3))
  148. sapply(c("png", "svg", "pdf"), function(fmt)
  149. ggsave(
  150. plot = reads_per_condition,
  151. filename = sprintf(
  152. "%s/%s/Supplementary Figure %i - Reads distribution per condition.%s",
  153. params$basic_outdir,
  154. params$microbiome_data_type,
  155. SuppFigIndex,
  156. fmt
  157. ),
  158. height = 8,
  159. width = 8,
  160. dpi = 600,
  161. bg = "white"
  162. )) %>% invisible()
  164. SuppFigIndex <- sum(SuppFigIndex, 1)
  165. ```
  167. ## Alpha diversity does not vary in relation to sequencing depth
  169. ```{r}
  170. alphaDivIndices <- c(
  171. "dbp_dominance",
  172. "observed_richness",
  173. "shannon_diversity",
  174. "gini"
  175. )
  177. MetaPhlAn.tse <- addAlpha(MetaPhlAn.tse, assay.type = params$assayTested, index = alphaDivIndices)
  178. ```
  180. ```{r}
  181. histograms_alpha <- lapply(alphaDivIndices,
  182. function(div) {
  183. stats_basic <- shapiro.test(as.data.frame(colData(MetaPhlAn.tse))[[div]])
  184. text <- paste0("Shapiro-Wilk test: W = ", round(stats_basic$statistic, 2), "; P-value = ", round(stats_basic$p.value, 3))
  185. plot <- ggplot(as.data.frame(colData(MetaPhlAn.tse)), aes(x = !!sym(div))) +
  186. geom_histogram() +
  187. labs(
  188. title = div,
  189. caption = text
  190. )
  192. }
  193. ) %>%
  194. ggarrange(plotlist = .)
  196. ```
  198. ```{r, tab.cap="Assessment of dependence of alpha diversity measures in relation to sequencing depth"}
  199. tmp <- sapply(alphaDivIndices, function(div) {
  200. lm(as.formula(paste(div, "~", "number_reads")), data = colData(MetaPhlAn.tse))
  202. }, USE.NAMES = TRUE, simplify = FALSE)
  204. table_results.df <- tmp %>% lapply(broom::tidy) %>%
  205. bind_rows(.id = "AlphaIndex") %>%
  206. filter(!grepl("Intercept", term)) %>%
  207. as_tibble()
  209. alpha_vs_reads_plot <- ggarrange(
  210. histograms_alpha,
  211. ggtexttable(table_results.df, rows = NULL),
  212. nrow = 2
  213. )
  215. sapply(c("png", "svg", "pdf"), function(fmt)
  216. ggsave(plot = alpha_vs_reads_plot, filename = sprintf(
  217. "%s/%s/Supplementary Figure %i - Alpha Indices vs number of reads.%s",
  218. params$basic_outdir,
  219. params$microbiome_data_type,
  220. SuppFigIndex,
  221. fmt
  222. ),
  223. height = 8,
  224. width = 8,
  225. dpi = 600,
  226. bg = "white"
  227. )
  228. ) %>% invisible()
  230. SuppFigIndex <- sum(SuppFigIndex, 1)
  231. ```
  233. ```{r, tab.cap="Assessment of dependence of alpha diversity measures in relation to sequencing depth"}
  234. table_results.df %>%
  235. kbl() %>%
  236. row_spec(which(table_results.df$p.value < 0.05), bold = TRUE) %>%
  237. kable_styling()
  238. ```
  240. # Main figure
  242. ## A - Stackplot relative abundances at Phylum Level
  244. ```{r}
  245. panelA_stackplot <- plotAbundance(
  246. x = MetaPhlAn.tse,
  247. assay.type = params$assayTested,
  248. group = "Phylum") +
  249. facet_grid(cols = vars(genotype_treatment_legend), scales = "free_x") +
  250. theme_light() +
  251. scale_x_discrete(labels = "x") +
  252. theme(
  253. panel.spacing.x = unit(0, "lines"),
  254. strip.background.x = element_rect(fill = "white", colour = "gray40", linewidth = unit(0.8, "mm")),
  255. panel.grid = element_blank(),
  256. panel.border = element_rect(colour = "gray40", linewidth = unit(0.8, units = "mm")),
  257. strip.text = element_text(color = "black"),
  258. axis.text.x = element_text(color = "white"),
  259. axis.ticks.x = element_blank(),
  260. strip.text.x = element_markdown()
  261. ) +
  262. labs(
  263. fill = "Phylum",
  264. x = " "
  265. )
  267. panelA_stackplot$data$colour_by <- gsub("p__", "", panelA_stackplot$data$colour_by)
  268. panelA_stackplot$data$colour_by <- gsub("_", " ", panelA_stackplot$data$colour_by)
  269. ```
  271. ## B - PCoA
  273. ```{r}
  274. set.seed(1234)
  275. dimRedName <- paste("Fp7", params$beta_dist, "mds", sep = "_")
  277. MetaPhlAn.tse <- mia::addMDS(
  278. MetaPhlAn.tse,
  279. assay.type = params$assayTested,
  280. method = params$beta_dist,
  281. ncomponents = 3, name = dimRedName
  282. )
  284. variances_explained <- attr(reducedDim(MetaPhlAn.tse, dimRedName), "eig")
  285. prop_variances_explained <- variances_explained/sum(variances_explained)
  287. Fp7_colData_for_PCoA.df <- cbind.data.frame(as.data.frame(colData(MetaPhlAn.tse)), reducedDim(MetaPhlAn.tse, dimRedName) %>% as.data.frame() %>% set_names(paste("MDS", 1:ncol(.))))
  289. panelB_PCoA <- Fp7_colData_for_PCoA.df %>%
  290. ggplot(aes(x = `MDS 1`, y = `MDS 2`, color = genotype_treatment_legend, fill = genotype_treatment_legend)) +
  291. geom_point(size = 2.5, stroke = 1, color = "black") +
  292. geom_point(size = 2.5) +
  293. ggConvexHull::geom_convexhull(show.legend = FALSE, alpha = 0.2) +
  294. scale_color_manual(values = overall_color_palette) +
  295. scale_fill_manual(values = overall_color_palette) +
  296. theme(legend.position = "top", legend.text = element_markdown()) +
  297. labs(color = "Genotype and Treatment",
  298. fill = "Genotype and Treatment",
  299. x = paste0("PCoA Axis 1 [", round(prop_variances_explained[1]*100, 1), "%]"),
  300. y = paste0("PCoA Axis 2 [", round(prop_variances_explained[2]*100, 2), "%]")
  301. )
  303. # panelB_PCoA
  304. ```
  306. ## C - Pairwise bray-curtis dissimilarities
  308. ### Statistical test (Bray-curtis)
  310. ```{r}
  311. dissimName <- paste("Fp7", params$beta_dist, sep = "_")
  312. MetaPhlAn.tse <- addDissimilarity(MetaPhlAn.tse, method = "bray", name = dissimName, assay.type = params$assayTested)
  313. mtx <- metadata(MetaPhlAn.tse)[[dissimName]]
  314. mtx[upper.tri(mtx,diag = TRUE)] <- NA
  316. pairwise_dissim_bray <- reshape2::melt(mtx) %>%
  317. filter(complete.cases(.)) %>%
  318. dplyr::inner_join(., as.data.frame(colData(MetaPhlAn.tse)) %>% transmute(Var1 = uuid, genotype_treatment_legend = genotype_treatment_legend), by = "Var1") %>%
  319. dplyr::inner_join(., as.data.frame(colData(MetaPhlAn.tse)) %>% transmute(Var2 = uuid, genotype_treatment_legend = genotype_treatment_legend), by = "Var2") %>%
  320. dplyr::select(Var1, Var2, contains(c("uuid", "genotype", "treatment")), value) %>%
  321. mutate(pairs = paste(genotype_treatment_legend.y, genotype_treatment_legend.x, sep = " vs ")) %>%
  322. filter(genotype_treatment_legend.x != genotype_treatment_legend.y)
  324. ```
  326. ```{r}
  327. set.seed(1234)
  328. library(vegan)
  329. pairwise_adonis_stats <- pairwise.adonis(x = as.dist(metadata(MetaPhlAn.tse)[[dissimName]]),
  330. factors = MetaPhlAn.tse@colData$genotype_treatment,
  331. perm = 9999,
  332. p.adjust.m = "BH"
  333. )
  335. ```
  337. ### Boxplot dissimilarities
  339. ```{r}
  340. pairwise_adonis_stats_asterisks <- pairwise_adonis_stats %>%
  341. mutate(pairs = gsub("Thy1-ASO F. prausnitzii", "Thy1-ASO *F. prausnitzii*", pairs))
  343. panelC_pairwise_dissim_boxplot <- pairwise_dissim_bray %>%
  344. ggplot(
  345. aes(x = pairs, y = value)
  346. ) +
  347. geom_point(position = position_jitter(width = 0.2, seed = 1234), stroke = 1.1) +
  348. geom_point(position = position_jitter(width = 0.2, seed = 1234), color = "gray30") +
  349. geom_boxplot(outliers = FALSE, alpha = 0.6, width = 0.3) +
  350. geom_text(data = pairwise_adonis_stats_asterisks %>%
  351. dplyr::left_join(
  352. group_by(pairwise_dissim_bray, pairs) %>%
  353. reframe(value = max(value)*1.1), by = "pairs") %>%
  354. mutate(pairs = factor(pairs, levels = pairs)),
  355. mapping = aes(x = pairs, y = value, label = sprintf("FDR = %.3f", p.adjusted)), hjust = 0, angle = 90) +
  356. # Apply the colorize_labels function to your x-axis labels
  357. scale_x_discrete(labels = function(x) colorize_labels(x, overall_color_palette)) +
  358. # Tell ggplot to render the axis text as Markdown/HTML
  359. theme(
  360. axis.text.x = ggtext::element_markdown(angle = 90, hjust = 1, vjust = 0.5)
  361. ) +
  362. ylim(min(pairwise_dissim_bray$value)*0.9, max(pairwise_dissim_bray$value)*1.5) +
  363. labs(
  364. x = "",
  365. y = "Bray-Curtis Pairwise Dissimilarity"
  366. )
  368. # make a horizontal version just in case
  370. panelC_pairwise_dissim_boxplot_horiz <- pairwise_dissim_bray %>%
  371. ggplot(
  372. aes(y = pairs, x = value)
  373. ) +
  374. geom_point(position = position_jitter(height = 0.2, seed = 1234), stroke = 1.1) +
  375. geom_point(position = position_jitter(height = 0.2, seed = 1234), color = "gray30") +
  376. geom_boxplot(outliers = FALSE, alpha = 0.6, width = 0.3) +
  377. geom_text(data = pairwise_adonis_stats_asterisks %>%
  378. dplyr::left_join(
  379. group_by(pairwise_dissim_bray, pairs) %>%
  380. reframe(value = max(value)*1.1), by = "pairs") %>%
  381. mutate(pairs = factor(pairs, levels = pairs)),
  382. mapping = aes(y = pairs, x = value, label = sprintf("FDR = %.3f", p.adjusted)), hjust = 0) +
  383. # Apply the colorize_labels function to your x-axis labels
  384. scale_y_discrete(labels = function(x) colorize_labels(x, overall_color_palette)) +
  385. # Tell ggplot to render the axis text as Markdown/HTML
  386. theme(axis.text.y = ggtext::element_markdown(angle = 0, hjust = 1)) +
  387. xlim(min(pairwise_dissim_bray$value)*0.9, max(pairwise_dissim_bray$value)*1.5) +
  388. labs(
  389. y = "",
  390. x = "Bray-Curtis Pairwise Dissimilarity"
  391. )
  392. ```
  394. ## D - Differentially abundant features due to treatment in Thy1-ASO mice
  396. Model formula is:
  398. $$
  399. \log_2(Microbiome)_{Thy1-ASO} \sim treament
  400. $$
  402. ```{r}
  403. # prepare input data
  404. tmp.tse <- MetaPhlAn.tse[,MetaPhlAn.tse@colData$genotype != "WT"]
  405. tmp.tse@colData$treatment <- relevel(tmp.tse@colData$treatment, "Control")
  406. tmp.tse@colData$treatment <- droplevels(tmp.tse@colData$treatment)
  407. tmp.tse <- tmp.tse[rowSums(assay(tmp.tse, params$assayTested)) > 0,]
  409. tmp.tse
  411. # define experiment name and output directories
  412. maaslin3_expmt <- "maaslin3_treatments_vs_Control_in_Thy1-ASO"
  413. maaslin3_final_outdir <- file.path(params$basic_outdir, params$microbiome_data_type, maaslin3_expmt)
  415. # create definitive output directory
  416. dir.create(maaslin3_final_outdir, showWarnings = FALSE, recursive = TRUE)
  417. ```
  419. Results are stored in ``r maaslin3_final_outdir``
  421. ```{r}
  422. set.seed(1234)
  423. treatments_vs_Control_in_Thy1ASO.maaslin3 <- maaslin3(
  424. input_data = tmp.tse,
  425. formula = ~ treatment,
  426. transform = "LOG",
  427. output = maaslin3_final_outdir,
  428. # do not standardize continuous metadata variables in the model
  429. standardize = FALSE,
  430. plot_associations = TRUE,
  431. save_plots_rds = TRUE,
  432. verbosity = "ERROR",
  433. summary_plot_first_n = 15
  434. )
  436. # enrich results with taxonomy information only in case of Metaphlan
  438. if(agrepl("metaphlan", params$assayTested)){
  439. ## add taxonomy to output tables: all_results.tsv
  440. dplyr::right_join(GTDB.df %>% mutate(SGB = paste0("t__", SGB)), read_tsv(file.path(maaslin3_final_outdir, "all_results.tsv")) %>% dplyr::rename(!!sym(params$microbiome_data_type) := feature), by = params$microbiome_data_type) %>%
  441. #arrange again based on increasing qval_individual
  442. arrange(qval_individual) %>%
  443. write_tsv(file.path(maaslin3_final_outdir, "all_results.tsv"))
  445. ## add taxonomy to output tables: significant_results.tsv
  446. dplyr::right_join(GTDB.df %>% mutate(SGB = paste0("t__", SGB)), read_tsv(file.path(maaslin3_final_outdir, "significant_results.tsv")) %>% dplyr::rename(!!sym(params$microbiome_data_type) := feature), by = params$microbiome_data_type) %>%
  447. write_tsv(file.path(maaslin3_final_outdir, "significant_results.tsv"))
  448. }
  449. # save raw output
  450. saveRDS(treatments_vs_Control_in_Thy1ASO.maaslin3, file = file.path(maaslin3_final_outdir, "maaslin3_raw_output.rds"))
  451. ```
  453. ```{r}
  454. GTDB.df <- GTDB.df %>%
  455. mutate(
  456. altNames = sapply(
  457. strsplit(paste(trimws(Genus), trimws(Species), sprintf("[%s]", SGB)), " ", fixed = TRUE),
  458. function(x) paste(unique(x), collapse = " ")
  459. )
  460. )
  462. altNames.chr <- GTDB.df$altNames
  463. names(altNames.chr) <- GTDB.df$SGB
  465. panelD_diffAbund <- readRDS(file = file.path(maaslin3_final_outdir, "figures", "summary_plot_gg.RDS")) +
  466. theme(strip.text = element_blank()) +
  467. # expand plot to show full errorbars
  468. scale_x_continuous(expand = expansion(mult = c(.1))) +
  469. # remove the null HP that would mess with the panel rendering
  470. guides(linetype = "none")
  474. panelD_diffAbund2 <- panelD_diffAbund +
  475. scale_y_discrete(
  476. labels = function(y) {
  477. altnamesBasic <- altNames.chr[gsub("t__", "", y)]
  478. altnamesItalics <- strsplit(altnamesBasic, "\\ ") %>% lapply(function(x)
  479. ifelse(str_ends(x, "cter|culum|terium|cter|monas|ccus|llus|cus|ctor|ella"), paste0("*", x, "*"), x)) %>%
  480. purrr::map_chr(.f = function(x) paste(x, collapse = " "))
  481. return(altnamesItalics)
  482. }
  483. ) + theme(
  484. axis.text.x = element_text(size = 12),
  485. axis.text.y = element_text(size = 12)
  486. )
  487. ```
  489. ## E - Mediation analysis with PCoA axes 1-3 on constipation and motor symptoms
  491. - Constipation symptoms: $Bead~expulsion~(minutes)$
  492. - Motor symptoms: $\frac{slips}{steps}~on~beam$
  494. We use the standard 3-equation approach:
  496. - (1) $Mediator \sim Treatment$
  498. - (2) $Outcome \sim Treatment$
  500. - (3) $Outcome \sim Mediator + Treatment$
  502. Where `Mediator` = PCo.x, `Outcome` = Bead Expulsion or Errors/Steps, `Treatment`
  503. *F. prausnitzii* supplementation.
  505. ```{r}
  506. # gather bead_expulsion_n data
  507. motor_gi_symptoms.list <- sapply(c("Beam_steps", "Bead_exp"), function(x) openxlsx::read.xlsx("Data/Fp7_GI_motor_assessment.xlsx", sheet = x), simplify = FALSE, USE.NAMES = TRUE)
  509. # prepare bead expulsion variable
  510. motor_gi_symptoms.list$Bead_exp <- motor_gi_symptoms.list$Bead_exp %>%
  511. dplyr::transmute(
  512. mouse_id = ID,
  513. bead_expulsion_min = Trial1)
  515. # prepare slips over steps variable
  516. motor_gi_symptoms.list$Beam_steps <- motor_gi_symptoms.list$Beam_steps %>%
  517. rowwise() %>%
  518. transmute(
  519. mouse_id = ID,
  520. Slips = Slips_Trial1,
  521. Steps = Steps_Trial1,
  522. # older choice of presenting data
  523. # Slips = mean(c(Slips_Trial1, Slips_Trial2), na.rm = TRUE),
  524. # Steps = mean(c(Steps_Trial1, Steps_Trial2), na.rm = TRUE),
  525. slips_over_steps = Slips/Steps)
  527. additional_variables.df <- purrr::reduce(motor_gi_symptoms.list, dplyr::full_join, by = "mouse_id")
  529. # incorporate these columns into the colData of the input
  530. new_df <- DataFrame(dplyr::left_join(as.data.frame(colData(MetaPhlAn.tse)), additional_variables.df, by = "mouse_id"))
  531. rownames(new_df) <- new_df$uuid
  533. colData(MetaPhlAn.tse) <- new_df
  534. ```
  536. ### Mediation analysis of bead expulsion
  538. ```{r}
  539. library(mediation)
  540. set.seed(1234)
  542. mediationData.df <- cbind.data.frame(
  543. as.data.frame(colData(MetaPhlAn.tse)),
  544. as.data.frame(reducedDim(MetaPhlAn.tse, dimRedName))) %>%
  545. dplyr::rename(
  546. PCo.1 = V1,
  547. PCo.2 = V2,
  548. PCo.3 = V3
  549. ) %>%
  550. filter(
  551. # remove NAs from covariates of interest
  552. complete.cases(dplyr::select(., contains("PCo"), genotype_treatment, bead_expulsion_min)),
  553. # remove WT Control and use ASO control as baseline
  554. genotype_treatment != "WT Control"
  555. ) %>%
  556. mutate(
  557. treatment = droplevels(treatment)
  558. )
  560. # make 3 mediation analyses as list, one per Principal Component
  561. mediations_Princ.Coord <- sapply(grepv("PCo", colnames(mediationData.df)), function(x) {
  562. # Model 1 - Effect of exposure on mediator
  563. mediation.lm <- lm(as.formula(paste(x, "treatment", sep = "~")), data = mediationData.df)
  564. # Model 2 - Effect of exposure on outcome (Total effect)
  565. outcome.lm <- lm(as.formula(paste("bead_expulsion_min ~ treatment", x, sep = "+")), data = mediationData.df)
  566. # Model 3 - Meidation analysis
  567. mediation_analysis = mediate(
  568. model.m = mediation.lm,
  569. model.y = outcome.lm,
  570. sims = 1000,
  571. treat = "treatment",
  572. mediator = x,
  573. boot = FALSE,
  574. robustSE = TRUE
  575. )
  576. # return all 3 models for completeness
  578. return(
  579. list(
  580. mediator.lm = mediation.lm,
  581. outcome.lm = outcome.lm,
  582. mediation_analysis = mediation_analysis
  583. )
  584. )
  585. }, simplify = FALSE, USE.NAMES=TRUE)
  587. data_as_long_table_bead_expulsion <- mediations_Princ.Coord %>%
  588. lapply("[[", "mediation_analysis") %>%
  589. lapply(tidy_mediation) %>%
  590. bind_rows(.id = "PCo") %>%
  591. mutate(Exposure = "Bead Expulsion (min)")
  593. ```
  595. ### Mediation analysis of Errors/steps
  597. ```{r}
  598. library(mediation)
  599. set.seed(1234)
  601. mediationData.df <- cbind.data.frame(
  602. as.data.frame(colData(MetaPhlAn.tse)),
  603. as.data.frame(reducedDim(MetaPhlAn.tse, dimRedName))) %>%
  604. dplyr::rename(
  605. PCo.1 = V1,
  606. PCo.2 = V2,
  607. PCo.3 = V3
  608. ) %>%
  609. filter(
  610. # remove NAs from covariates of interest
  611. complete.cases(dplyr::select(., contains("PCo"), genotype_treatment, slips_over_steps)),
  612. # remove WT Control and use ASO control as baseline
  613. genotype_treatment != "WT Control"
  614. ) %>%
  615. mutate(
  616. treatment = droplevels(treatment)
  617. )
  619. # make 3 mediation analyses as list, one per Principal Component
  620. mediations_Princ.Coord <- sapply(grepv("PCo", colnames(mediationData.df)), function(x) {
  621. # Model 1 - Effect of exposure on mediator
  622. mediation.lm <- lm(as.formula(paste(x, "treatment", sep = "~")), data = mediationData.df)
  623. # Model 2 - Effect of exposure on outcome (Total effect)
  624. outcome.lm <- lm(as.formula(paste("slips_over_steps ~ treatment", x, sep = "+")), data = mediationData.df)
  625. # Model 3 - Meidation analysis
  626. mediation_analysis = mediate(
  627. model.m = mediation.lm,
  628. model.y = outcome.lm,
  629. sims = 1000,
  630. treat = "treatment",
  631. mediator = x,
  632. boot = FALSE,
  633. robustSE = TRUE
  634. )
  635. # return all 3 models for completeness
  637. return(
  638. list(
  639. mediator.lm = mediation.lm,
  640. outcome.lm = outcome.lm,
  641. mediation_analysis = mediation_analysis
  642. )
  643. )
  644. }, simplify = FALSE, USE.NAMES=TRUE)
  646. data_as_long_table_slips_over_steps <- mediations_Princ.Coord %>%
  647. lapply("[[", "mediation_analysis") %>%
  648. lapply(tidy_mediation) %>%
  649. bind_rows(.id = "PCo") %>%
  650. mutate(Exposure = "Beam Traversal - Errors/Steps")
  651. ```
  653. ### Panel E plot generation
  655. ```{r}
  656. data_as_long_table_for_horizPlot <- rbind.data.frame(data_as_long_table_slips_over_steps,data_as_long_table_bead_expulsion) %>%
  657. mutate(
  658. PCo = factor(PCo, levels = sort(unique(data_as_long_table_slips_over_steps$PCo), decreasing = TRUE)),
  659. Effect = factor(Effect, levels = sort(unique(Effect), decreasing = TRUE))
  660. )
  661. ```
  663. ```{r}
  664. panelE_mediation_motor_GI_symptoms <- ggplot(data_as_long_table_for_horizPlot, aes(x = Estimate, y = Effect, color = PCo, shape = P.value < 0.05)) +
  665. geom_vline(xintercept = 0, lty = "dashed", color = "darkgray") +
  666. geom_errorbar(aes(xmin = CI_Lower_95, xmax = CI_Upper_95, y = Effect, color = PCo), position = position_dodge(width = 0.5, preserve = "total"), width = 0, show.legend = FALSE) +
  667. geom_point(position = position_dodge(width = 0.5, preserve = "total"), size = 4) +
  668. scale_shape_manual(values = c(1,19)) +
  669. facet_nested(cols = vars(Exposure), scales = "free_x") +
  670. scale_color_manual(values = c("darksalmon", "brown3", "brown4"), labels = c("PCo 3", "PCo 2", "PCo 1")) +
  671. guides(color = guide_legend(reverse = TRUE), shape = "none") +
  672. theme(
  673. panel.spacing.x = unit(0, "lines"),
  674. strip.background.x = element_rect(fill = "white", colour = "gray40", linewidth = unit(0.8, "mm")),
  675. strip.text = element_text(color = "black"),
  676. axis.text = element_text(size = 12),
  677. axis.title = element_text(size = 12)
  678. ) +
  679. labs(
  680. color = NULL,
  681. y = NULL,
  682. x = expression(beta~Estimate)
  683. )
  685. panelE_mediation_motor_GI_symptoms
  687. panelE_mediation_motor_GI_symptoms_shorterLabels <- panelE_mediation_motor_GI_symptoms + scale_y_discrete(labels = function(x) gsub("Average\\ |Causal\\ ", "", x))
  690. ```
  692. ## Arrange panels and save figure
  694. ```{r, fig.height=8, fig.width=12}
  696. panel_A_legend <- get_legend(panelA_stackplot + guides(fill = guide_legend(nrow = 2, title = "Phylum", position = "top")))
  698. Microbiome_Main_Figure <- ggarrange(
  700. ggarrange(
  701. panelA_stackplot + guides(fill = guide_legend(title = "Phylum")) + labs(y = "Relative Abundance", x = "Mouse") + theme(axis.title = element_text(size = 12)),
  702. panelB_PCoA + theme(axis.title = element_text(size = 12)),
  703. common.legend = TRUE, widths = c(1, 1), legend = "top", labels = c("A", "B")),
  705. ggarrange(
  706. panelC_pairwise_dissim_boxplot_horiz + theme(axis.title = element_text(size = 12), axis.text = element_text(size = 9)),
  707. ggarrange(
  708. get_legend(panelB_PCoA + guides(color = guide_legend(ncol = 1, position = "right"))),
  709. ggplot() + theme_void(), nrow = 2, heights = c(1, 0.35)
  710. ),
  712. nrow = 1, widths = c(1,0.2), labels = c("C", "")
  713. ),
  715. ggarrange(
  716. ggarrange(
  717. ggplot() + theme_void(),
  718. panelD_diffAbund2 +
  719. coord_flip() +
  720. theme(axis.text.x = element_markdown(angle = 90, hjust = 1, vjust = 0.5), axis.title = element_text(hjust = 1),
  721. legend.position = "none") +
  722. labs(x = expression(log[2]~Fold~Change), y = NULL) +
  723. theme(axis.title = element_text(size = 12)),
  724. panelE_mediation_motor_GI_symptoms_shorterLabels + theme(legend.position = "bottom"),
  725. widths = c(0.05, 0.85, 1), nrow = 1, labels = c("D", "", "E")
  726. ),
  727. get_legend(panelD_diffAbund2 + theme(legend.position = "bottom", legend.text = element_text(size = 10), legend.title = element_text(size = 12))), nrow = 2,
  728. heights = c(1, 0.2)
  729. ),
  731. nrow = 3,
  732. heights = c(1, 0.25, 1)
  733. )
  735. sapply(c("png", "svg", "pdf"), function(fmt)
  736. ggsave(plot = Microbiome_Main_Figure, filename = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Main_Microbiome_Figure.%s", fmt)), width = 11, height = 13, bg = "white")
  737. )
  738. ```
  740. # Supplementary Figure `r SuppFigIndex` - alpha and F/B ratios vs genotype treatment
  742. ```{r}
  743. alpha_data.df <- MetaPhlAn.tse %>%
  744. colData() %>%
  745. as.data.frame() %>%
  746. dplyr::rename(
  747. `Shannon Diversity` = shannon_diversity,
  748. `Observed Richness` = observed_richness
  749. )
  751. comparisons <- combn(levels(alpha_data.df$genotype_treatment_legend), 2, simplify = FALSE)
  752. ```
  754. ## A - Observed Richness
  756. ```{r}
  757. anova_test <- broom::tidy(aov(`Observed Richness` ~ genotype_treatment, data = alpha_data.df)) %>%
  758. as.data.frame()
  760. anova_test.text <- sprintf("ANOVA: df = %i, F = %.3f, P = %.4f", anova_test[1,"df"],
  761. anova_test[1,"statistic"], anova_test[1,"p.value"])
  763. panelA_Obs_boxplot <- alpha_data.df %>%
  764. ggplot(aes(x = genotype_treatment_legend, y = `Observed Richness`, color = genotype_treatment_legend, fill = genotype_treatment_legend)) +
  765. geom_boxplot(outliers = FALSE, show.legend = FALSE, color = "black", alpha = 0.4) +
  766. geom_point(color = "black", stroke = 1.1, position= position_jitterdodge(jitter.width = 0.5, seed = 1234)) +
  767. geom_point(position= position_jitterdodge(jitter.width = 0.5, seed = 1234)) +
  768. scale_color_manual(values = overall_color_palette) +
  769. scale_fill_manual(values = overall_color_palette) +
  770. # stat_pvalue_manual(stats_obs, label = "p.adj.round3", tip.length = 0.01) +
  771. theme(axis.text.x = element_blank(), legend.text = element_markdown()) +
  772. labs(
  773. x = NULL,
  774. color = "Genotype and Treatment",
  775. fill = "Genotype and Treatment",
  776. caption = anova_test.text)
  777. ```
  779. ## B - Shannon Diversity
  781. ```{r}
  782. anova_test <- broom::tidy(aov(`Shannon Diversity` ~ genotype_treatment, data = alpha_data.df)) %>%
  783. as.data.frame()
  785. anova_test.text <- sprintf("ANOVA: df = %i, F = %.3f, P = %.4f", anova_test[1,"df"],
  786. anova_test[1,"statistic"], anova_test[1,"p.value"])
  788. panelB_Shannon_boxplot <- alpha_data.df %>%
  789. ggplot(aes(x = genotype_treatment_legend, y = `Shannon Diversity`, color = genotype_treatment_legend, fill = genotype_treatment_legend)) +
  790. geom_boxplot(outliers = FALSE, show.legend = FALSE, color = "black", alpha = 0.4) +
  791. geom_point(color = "black", stroke = 1.1, position= position_jitterdodge(jitter.width = 0.5, seed = 1234)) +
  792. geom_point(position= position_jitterdodge(jitter.width = 0.5, seed = 1234)) +
  793. scale_color_manual(values = overall_color_palette) +
  794. scale_fill_manual(values = overall_color_palette) +
  795. # stat_pvalue_manual(stats_shannon, label = "p.adj.round3", tip.length = 0.01) +
  796. theme(axis.text.x = element_blank(), legend.text = element_markdown()) +
  797. labs(
  798. x = NULL,
  799. color = "Genotype and Treatment",
  800. fill = "Genotype and Treatment",
  801. caption = anova_test.text
  802. )
  803. ```
  805. ## C - Firmicutes/Bacteroidetes Ratio
  807. ```{r}
  808. assayextractedPhylum <- assay(agglomerateByRank(MetaPhlAn.tse, "Phylum"), params$assayTested)
  809. MetaPhlAn.tse$FB.Ratio <- assayextractedPhylum["p__Firmicutes",]/assayextractedPhylum["p__Bacteroidota",]
  811. anova_test <- broom::tidy(aov(FB.Ratio ~ genotype_treatment, data = colData(MetaPhlAn.tse))) %>%
  812. as.data.frame()
  813. anova_test.text <- sprintf("ANOVA: df = %i, F = %.3f, P = %.4f", anova_test[1,"df"],
  814. anova_test[1,"statistic"], anova_test[1,"p.value"])
  816. panelF_FB.Ratio.df <- colData(MetaPhlAn.tse) %>%
  817. as.data.frame() %>%
  818. dplyr::select(genotype_treatment, genotype_treatment_legend, FB.Ratio)
  820. panelC_FB.Ratio <- panelF_FB.Ratio.df %>%
  821. ggplot(aes(x = genotype_treatment_legend, y = FB.Ratio, color = genotype_treatment_legend, fill = genotype_treatment_legend)) +
  822. geom_boxplot(outliers = FALSE, show.legend = FALSE, color = "black", alpha = 0.4) +
  823. geom_point(color = "black", stroke = 1.1, position= position_jitterdodge(jitter.width = 0.5, seed = 1234)) +
  824. geom_point(position= position_jitterdodge(jitter.width = 0.5, seed = 1234)) +
  825. scale_color_manual(values = overall_color_palette) +
  826. scale_fill_manual(values = overall_color_palette) +
  827. # stat_pvalue_manual(stats_fbrartio, label = "p.adj.round3", tip.length = 0.01) +
  828. theme(axis.text.x = element_blank(), legend.text = element_markdown()) +
  829. labs(
  830. x = NULL,
  831. y = "Firmicutes/Bacteroidetes Ratio",
  832. color = "Genotype and Treatment",
  833. fill = "Genotype and Treatment",
  834. caption = anova_test.text)
  836. ```
  838. ## Arrange panels and save figure
  840. ```{r}
  841. microbiome_panel_Supplementary <- ggarrange(
  842. panelA_Obs_boxplot + theme(axis.ticks.x = element_blank()) + scale_y_continuous(expand = expansion(mult = 0.08)),
  843. panelB_Shannon_boxplot + theme(axis.ticks.x = element_blank()) + scale_y_continuous(expand = expansion(mult = 0.08)),
  844. panelC_FB.Ratio + theme(axis.ticks.x = element_blank()) + scale_y_continuous(expand = expansion(mult = 0.08)),
  845. nrow = 1, labels = c("A", "B", "C"),
  846. common.legend = TRUE)
  848. sapply(c("png", "svg", "pdf"), function(fmt)
  849. ggsave(plot = microbiome_panel_Supplementary, filename = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Figure %i - Alpha richness, diversity and FB Ratio.%s", SuppFigIndex, fmt)), width = 10, height = 4)
  850. )
  851. SuppFigIndex <- sum(SuppFigIndex, 1)
  853. microbiome_panel_Supplementary
  854. ```
  856. # Supplementary Tables
  858. ## Supplementary Table `r SuppTabIndex` - BugSigDB Parkinson gut microbiome siganture across studies
  860. ```{r}
  861. host <- "Homo sapiens"
  863. bsdb <- importBugSigDB(version = "devel", cache=FALSE)
  865. bsdb_PD_signatures_allMethods <- filter(bsdb,
  866. # get studies/signatures on Parkinson alone
  867. grepl("Parkinson's disease", Condition),
  868. # restrict search to Homo sapiens
  869. `Host species` == "Homo sapiens"
  870. ) %>%
  871. # fix taxa names in case k__ is duplicated because of some weird
  872. # changes found on May 20 2025
  873. mutate(
  874. `MetaPhlAn taxon names` = lapply(`MetaPhlAn taxon names`, function(x) gsub("|k__", "_", x = x, fixed = TRUE))
  875. )
  877. # fix typo of feces,blood
  879. bsdb_PD_signatures_allMethods <- mutate(bsdb_PD_signatures_allMethods,
  880. `Body site` = ifelse(`Body site` == "Feces,Blood", "Blood,Feces", `Body site`))
  882. bsdb_PD_signatures_allMethods_onlyFeces <- filter(bsdb_PD_signatures_allMethods, `Body site` == "Feces")
  884. top50_bugsigdb_PD_taxa <- createTaxonTable(bsdb_PD_signatures_allMethods_onlyFeces, n = 50) %>%
  885. filter(`Taxonomic Level` %in% c("genus", "species")) %>%
  886. arrange(`Binomial Test pval`) %>%
  887. dplyr::select(`Taxon Name`:`Binomial Test pval`, metaphlan_name)
  889. write_tsv(top50_bugsigdb_PD_taxa, file = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Table %i - BugSigDB Parkinson disease signature as of %s.tsv", SuppTabIndex, Sys.Date())))
  891. SuppTabIndex <- SuppTabIndex + 1
  892. ```
  894. ## Supplementary Table `r SuppTabIndex` - PERMANOVA relative to Main Figure Panel B
  896. ```{r}
  897. MetaPhlAn.tse <- addPERMANOVA(MetaPhlAn.tse,
  898. formula = x ~ genotype_treatment,
  899. assay.type = params$assayTested,
  900. method = params$beta_dist,
  901. permutations = 9999,
  902. by = "margin")
  904. permanova_table.df <- PERMANOVA_to_table(MetaPhlAn.tse)
  906. write_tsv(permanova_table.df, file = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Table %i - PERMANOVA of %s distance relative to Main Figure panel B.tsv", SuppTabIndex, params$beta_dist)))
  908. SuppTabIndex <- SuppTabIndex + 1
  909. ```
  911. ## Supplementary Table `r SuppTabIndex` - Pairwise PERMANOVA statistics related to Main Figure Panel C
  913. ```{r}
  914. write_tsv(pairwise_adonis_stats, file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Table %i - pairwise PERMANOVA of %s distance relative to Main Figure panel C.tsv", SuppTabIndex, params$beta_dist)))
  916. SuppTabIndex <- SuppTabIndex + 1
  917. ```
  919. ## Supplementary Table `r SuppTabIndex` - Genome Taxonomy Database of the SGBs in the study
  921. ```{r}
  922. taxaNames_no_UN_no_Prefix <- gsub("t__", "", rownames(MetaPhlAn.tse)) %>%
  923. .[!grepl("^UN", .)]
  925. GTDB.df[match(taxaNames_no_UN_no_Prefix, GTDB.df$SGB),] %>%
  926. write_tsv(file = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Table %i - GTDB taxonomy of the SGBs in the study.tsv", SuppTabIndex)))
  928. SuppTabIndex <- SuppTabIndex + 1
  929. ```
  931. ## Supplementary Table `r SuppTabIndex` - Full Differential abundance analysis results
  933. ```{r}
  934. file.copy(from = file.path(params$basic_outdir, params$microbiome_data_type, maaslin3_expmt, "all_results.tsv"), to = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Table %i - MaAsLin3 Differential abundance analysis ASO_Treated vs ASO_Control.tsv", SuppTabIndex)))
  936. SuppTabIndex <- sum(SuppTabIndex, 1)
  937. ```
  939. ## Supplementary Table `r SuppTabIndex` - Mediation analysis on bead expulsion related to Main Figure panel E
  941. ```{r}
  942. write_tsv(data_as_long_table_bead_expulsion, file = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Table %i - Mediation analysis on bead expulsion against PCo 1,2,3.tsv", SuppTabIndex)))
  944. SuppTabIndex <- sum(SuppTabIndex, 1)
  945. ```
  947. ## Supplementary Table `r SuppTabIndex` - Mediation analysis on slips.over.steps ratio related to Main Figure panel E
  949. ```{r}
  950. write_tsv(data_as_long_table_bead_expulsion, file = file.path(params$basic_outdir, params$microbiome_data_type, sprintf("Supplementary Table %i - Mediation analysis on slips.over.steps ratio against PCo 1,2,3.tsv", SuppTabIndex)))
  952. SuppTabIndex <- sum(SuppTabIndex, 1)
  953. ```
  955. # Session Info
  957. ```{r}
  958. sessionInfo()
  959. ```

Main-Analysis-SGB.Rmd, no license · at the source

Overview

Authors: Anastasiya Moiseyenko1,2, Giacomo Antonello2,3,4, Aubrey M Schonhoff1,2, Joseph C Boktor1,2, Kaelyn Long2,3,4, Blake Dirks5, Anastasiya D Oguienko1, Alexander Viloria Winnett1,6, Patrick Simpson1, Dorsa Daeizadeh5, Rustem F Ismagilov1,7, Rosa Krajmalnik-Brown5,8, Nicola Segata2,4, Levi D Waldron2,3,4, Sarkis K Mazmanian1,2
ORCID iDs: Giacomo Antonello, Kaelyn Long, Dorsa Daeizadeh
  1. Division of Biology and Biological Engineering, California Institute of Technology, Pasadena, CA USA
  2. Aligning Science Across Parkinson’s (ASAP) Collaborative Research Network, Chevy Chase, MD USA
  3. Department of Epidemiology and Biostatistics and Institute for Implementation Science in Population Health, Graduate School of Public Health and Health Policy, City University of New York (CUNY), New York, NY USA
  4. Department of Cellular, Computational and Integrative Biology (CIBIO), University of Trento, Povo, Trento Italy
  5. Biodesign Center for Health through Microbiomes, Arizona State University, Tempe, AZ USA
  6. UCLA-Caltech Medical Scientist Training Program, University of California Los Angeles, Los Angeles, CA USA
  7. Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA USA
  8. School of Sustainable Engineering and the Built Environment, Arizona State University, Tempe, AZ USA
Institutions: California Institute of Technology (United States); Aligning Science Across Parkinson's (United States); City University of New York (United States); University of Trento (Italy); CUNY Graduate School of Public Health and Health Policy (United States); Arizona State University (United States); University of California, Los Angeles (United States)
Journal: NPJ Parkinson's disease, volume 12, issue 1, article 94
Dates: received 18 September 2025; accepted 30 January 2026; published online 5 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41531-026-01287-x · PMID 41786774 · PMCID PMC13077053 · OpenAlex W7133788201
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), mouse (organism), Parkinson's (population), cellular / molecular (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Physiology & signal measures
Keywords: Microbiology, Neurology, Neuroscience
Topic: Parkinson's Disease Mechanisms and Treatments (Neurology, Medicine), according to OpenAlex
Funding: Aligning Science Across Parkinson's (ASAP-020495, ASAP-000375); Heritage Medical Research Institute; Aligning Science Across Parkinson&apos;s (ASAP-020495, ASAP-000375); Beckman Institute at Caltech; NCI NIH HHS (R01 CA230551); Jacobs Institute for Molecular Engineering for Medicine
Citations: cited by 2 papers (Europe PMC); 152 references in the paper
Research resources: RRID:AB_2099233, RRID:AB_2714215, 1:1000) or Anti-β-actin Mouse RRID:AB_306371, RRID:AB_330924, RRID:AB_763537, RRID:AB_869973, in R RRID:SCR_001905, RRID:SCR_002798, RRID:SCR_004463, RRID:SCR_004915, RRID:SCR_008394, RRID:SCR_008520, RRID:SCR_012919, was performed using ggplot2 RRID:SCR_014601, HUMANnN v3.9 RRID:SCR_014620, RRID:SCR_015687, RRID:SCR_016884, B and 2A were created in BioRender RRID:SCR_018361, RRID:SCR_023241, The mia package RRID:SCR_023619, RRID:SCR_026984

Abstract

Gut microbiome composition is altered in Parkinson’s disease (PD), the fastest-growing neurological condition, that is characterized by neurodegeneration, motor dysfunction, and is frequently accompanied by gastrointestinal (GI) symptoms. Notably, microbial taxa with anti-inflammatory properties are consistently depleted in PD patients compared to controls. To explore whether specific gut bacteria may be disease-protective, we assembled a microbial consortium of 8 human-associated taxa that are reduced in individuals with PD. Treatment of α-synuclein overexpressing (Thy1-ASO) mice, an animal model of PD, with this consortium improved motor and GI deficits. A single bacterial species from this consortium, Faecalibacterium prausnitzii, was sufficient to correct gut microbiome deviations in Thy1-ASO mice, induce anti-inflammatory immune responses, and promote protective colonic gene expression profiles. Accordingly, oral treatment with F. prausnitzii robustly ameliorated motor and GI deficits and reduced α-synuclein aggregates in the brain. These findings support the emerging hypothesis of functional contributions by the microbiome to PD outcomes, and embolden the development of potential probiotic therapies to treat motor and non-motor symptoms.

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

Repositories

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

pmartinezarbizu/pairwiseAdonis

License: GPL-3.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: cb190f7668a0c82c0b0853927db239e7b9ec3e83, 2 June 2023
Languages: R (2)
Size: 12 files, 2 scripts
Software Heritage: not archived
Found in: the text, “Fecal microbiome analysis with shotgun metagenom”
Holds: README, license file, environment (pairwiseAdonis/DESCRIPTION), documentation
Not found: CITATION.cff, tests, continuous integration
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
4 files
At the source: github.com/pmartinezarbizu/pairwiseAdonis

Zenodo 17137678

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
At the source:

Zenodo 17128213

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Languages: R (3)
Size: 5 files, 3 scripts
Software Heritage: not checked
Found in: “Data availability”
Holds: README, 3 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (3 files), ggpubr (2 files), broom (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
4 files
At the source:

Zenodo 17981510

License: CC-BY-4.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: clusterProfiler (2 files), DESeq2 (2 files), tidyverse (2 files), FastQC (1 file), ggpubr (1 file), patchwork (1 file), pheatmap (1 file), SAMtools (1 file), STAR (1 file), Subread (featureCounts) (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
4 files
At the source:

anamois/amoiseyenko2025

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 521ea8c08870f321f01eca853104c64930cbe4b4, 14 February 2026
Languages: Jupyter (4)
Size: 289 files, 4 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, license file, 4 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (4 files), pandas (4 files), SciPy (4 files), scikit-posthocs (2 files)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
6 files
At the source: github.com/anamois/amoiseyenko2025/tree/v4

jboktor/fprausnitzii-treatment-rnaseq

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 66d3461a83022c15e6dc032bad5b8148ba0d2cf3, 18 December 2025
Languages: R (1), Quarto (1), Shell (1)
Size: 8 files, 3 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, 1 notebook
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: clusterProfiler (2 files), DESeq2 (2 files), tidyverse (2 files), FastQC (1 file), ggpubr (1 file), patchwork (1 file), pheatmap (1 file), SAMtools (1 file), STAR (1 file), Subread (featureCounts) (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
4 files
At the source: github.com/jboktor/fprausnitzii-treatment-rnaseq/tree/v1.0.1

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:

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

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

Data

Datasets cited

Data availability

The data, code, protocols, and key lab materials used and generated in this study are listed in the Key Resource Table (Supplementary Table 2, also available at 10.5281/zenodo.17137678) alongside their persistent identifiers. For microbiome analysis data, raw FASTQ files were deposited in the SRA (PRJNA1259538) along with sample metadata, and code to generate data, figures, and statistics is available at 10.5281/zenodo.17128213. For RNA sequencing analysis data, raw FASTQ files were deposited in the SRA (PRJNA1308739) along with sample metadata, and code to generate data, figures, and statistics is available at 10.5281/zenodo.17981510. Raw flow cytometry data were deposited to Zenodo (10.5281/zenodo.16929959) along with sample metadata. All other raw data and analyses, including behavior data analysis, original blot images, multiplex, and VFA data, can be found at 10.5281/zenodo.17137678.

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, 15 authors, 3 keywords, 6 funders, 144 references, 21 RRIDs.

Cite

This paper

Moiseyenko, A., Antonello, G., Schonhoff, A. M., Boktor, J. C., Long, K., Dirks, B., Oguienko, A. D., Winnett, A. V., Simpson, P., Daeizadeh, D., Ismagilov, R. F., Krajmalnik-Brown, R., Segata, N., Waldron, L. D., & Mazmanian, S. K. (2026). Faecalibacterium prausnitzii, depleted in the Parkinson's disease microbiome, improves motor deficits in α-synuclein overexpressing mice. NPJ Parkinson's disease, 12(1), 94. https://doi.org/10.1038/s41531-026-01287-x

BibTeX

@article{moiseyenko2026faecalibacterium,
author = {Moiseyenko, Anastasiya and Antonello, Giacomo and Schonhoff, Aubrey M and Boktor, Joseph C and Long, Kaelyn and Dirks, Blake and Oguienko, Anastasiya D and Winnett, Alexander Viloria and Simpson, Patrick and Daeizadeh, Dorsa and Ismagilov, Rustem F and Krajmalnik-Brown, Rosa and Segata, Nicola and Waldron, Levi D and Mazmanian, Sarkis K},
title = {{Faecalibacterium prausnitzii, depleted in the Parkinson's disease microbiome, improves motor deficits in α-synuclein overexpressing mice}},
journal = {NPJ Parkinson's disease},
year = {2026},
month = mar,
volume = {12},
number = {1},
pages = {94},
publisher = {Nature Publishing Group},
issn = {2373-8057},
doi = {10.1038/s41531-026-01287-x},
url = {https://doi.org/10.1038/s41531-026-01287-x},
pmid = {41786774},
pmcid = {PMC13077053}
}

RIS

TY - JOUR
AU - Moiseyenko, Anastasiya
AU - Antonello, Giacomo
AU - Schonhoff, Aubrey M
AU - Boktor, Joseph C
AU - Long, Kaelyn
AU - Dirks, Blake
AU - Oguienko, Anastasiya D
AU - Winnett, Alexander Viloria
AU - Simpson, Patrick
AU - Daeizadeh, Dorsa
AU - Ismagilov, Rustem F
AU - Krajmalnik-Brown, Rosa
AU - Segata, Nicola
AU - Waldron, Levi D
AU - Mazmanian, Sarkis K
TI - Faecalibacterium prausnitzii, depleted in the Parkinson's disease microbiome, improves motor deficits in α-synuclein overexpressing mice
T2 - NPJ Parkinson's disease
J2 - NPJ Parkinsons Dis
PY - 2026
DA - 2026/03/05
VL - 12
IS - 1
SP - 94
SN - 2373-8057
PB - Nature Publishing Group
DO - 10.1038/s41531-026-01287-x
UR - https://doi.org/10.1038/s41531-026-01287-x
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41531-026-01287-x",
"type": "article-journal",
"title": "Faecalibacterium prausnitzii, depleted in the Parkinson's disease microbiome, improves motor deficits in α-synuclein overexpressing mice",
"container-title": "NPJ Parkinson's disease",
"author": [
{
"family": "Moiseyenko",
"given": "Anastasiya"
},
{
"family": "Antonello",
"given": "Giacomo"
},
{
"family": "Schonhoff",
"given": "Aubrey M"
},
{
"family": "Boktor",
"given": "Joseph C"
},
{
"family": "Long",
"given": "Kaelyn"
},
{
"family": "Dirks",
"given": "Blake"
},
{
"family": "Oguienko",
"given": "Anastasiya D"
},
{
"family": "Winnett",
"given": "Alexander Viloria"
},
{
"family": "Simpson",
"given": "Patrick"
},
{
"family": "Daeizadeh",
"given": "Dorsa"
},
{
"family": "Ismagilov",
"given": "Rustem F"
},
{
"family": "Krajmalnik-Brown",
"given": "Rosa"
},
{
"family": "Segata",
"given": "Nicola"
},
{
"family": "Waldron",
"given": "Levi D"
},
{
"family": "Mazmanian",
"given": "Sarkis K"
}
],
"container-title-short": "NPJ Parkinsons Dis",
"volume": "12",
"issue": "1",
"page": "94",
"DOI": "10.1038/s41531-026-01287-x",
"PMID": "41786774",
"PMCID": "PMC13077053",
"ISSN": "2373-8057",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41531-026-01287-x",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
5
]
]
}
}

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.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: FastQC, Subread (featureCounts), STAR, 12 other tools, 4 references
[2] doi:10.21203/rs.3.rs-9927928/v1 [code]
Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: FastQC, Subread (featureCounts), STAR, 11 other tools, mouse, 3 references
[3] doi:10.1126/sciadv.aed2952 [code]
Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.
Journal: Science advances
In common: Subread (featureCounts), STAR, SAMtools, 10 other tools, Parkinson's, cellular / molecular, 3 references
[4] doi:10.1038/s41586-026-10512-9 [code]
Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.
Journal: Nature
In common: Subread (featureCounts), STAR, SAMtools, 9 other tools, mouse, cellular / molecular, 3 references
[5] doi:10.1038/s41386-026-02406-1 [code]
Functional genomic profiling of schizophrenia-associated genes reveals key microglial regulators.
Journal: Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology
In common: FastQC, STAR, SAMtools, 7 other tools, cellular / molecular, 3 references
[6] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: Subread (featureCounts), SAMtools, broom, 10 other tools, mouse, cellular / molecular, 1 reference
[7] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Subread (featureCounts), SAMtools, broom, 10 other tools, mouse, cellular / molecular
[8] doi:10.7554/elife.107393 [code]
Chromosome-scale genome assembly of the European common cuttlefish &lt;i&gt;Sepia officinalis&lt;/i&gt;.
Journal: eLife
In common: Subread (featureCounts), STAR, SAMtools, 7 other tools, cellular / molecular, 4 references
[9] doi:10.1038/s41467-026-71432-w [code]
MeCP2 gene dosage-dependent neurodevelopmentally restricted defects arise by aberrant activation of cell fate-determining bivalent genes.
Journal: Nature communications
In common: FastQC, Subread (featureCounts), STAR, 8 other tools, mouse, cellular / molecular
[10] doi:10.1038/s41380-026-03578-4 [code]
Assessing molecular gene by treatment interactions using a population of neural progenitors exposed to valproic acid and lithium.
Journal: Molecular psychiatry
In common: FastQC, Subread (featureCounts), STAR, 6 other tools, cellular / molecular, 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.