From c5308a5072f0ca5771417ea7c2e20679728a68db Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Tue, 15 Jul 2025 01:36:24 +0800 Subject: [PATCH 01/16] Create multimedia.qmd --- inst/pages/multimedia.qmd | 590 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 590 insertions(+) create mode 100644 inst/pages/multimedia.qmd diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd new file mode 100644 index 00000000..439b2666 --- /dev/null +++ b/inst/pages/multimedia.qmd @@ -0,0 +1,590 @@ +--- +title: "Multimodal Mediation Analysis" +format: html +editor: visual +author: Yihan Liu, Himel Mallick +--- + +# **Multimodal Mediation Analysis** {#sec-MSEA} + +```{r setup, echo=FALSE, results="asis"} +remove(list = ls()) +invisible(gc()) +library(rebook) +chapterPreamble() + +``` + +Building upon the concepts and workflows presented in the previous chapter on +mediation analysis, here we demonstrate how to perform **multimodal mediation** +**analysis** to identify potential mediators across multiple omics layers. +**Multimodal mediation analysis** examines whether and how an exposure (e.g., +treatment, diet) affects an outcome (e.g., disease, phenotype) **through** +**intermediate variables (mediators) across multiple omics layers** +**simultaneously** (e.g., microbiome, metabolomics). It identifies which +features across modalities act as mediators, quantifies indirect (mediated) +effects while adjusting for covariates, and compares mediation strength across +layers, providing mechanistic insights into **how exposures impact outcomes** +**via molecular pathways**. + +```{r} +#| label: fig-multimodal-mediation +#| fig-cap: Directed acyclic graph illustrating multimodal mediation where an exposure affects multiple mediators across modalities (microbial species and metabolites), which in turn affect the outcome. A direct path from exposure to outcome is also included (dashed). While gut microbial species biologically influence metabolite levels, for tractability in high-dimensional mediation analysis we assume species and metabolites act as conditionally independent parallel mediators given exposure, acknowledging this as an approximation +#| echo: false + +library(DiagrammeR) + +grViz(" +digraph multimodal_mediation { + graph [layout = dot, rankdir = LR] + + node [fontname = Helvetica, fontsize = 12] + + Exposure [shape=box, style=filled, fillcolor=lightblue, color=black] + Outcome [shape=box, style=filled, fillcolor=palegreen, color=black] + + M1 [label='Species 1', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] + M2 [label='Species 2', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] + M3 [label='Metabolite 1', shape=ellipse, style=filled, fillcolor=lightpink, color=black] + M4 [label='Metabolite 2', shape=ellipse, style=filled, fillcolor=lightpink, color=black] + + Exposure -> M1 + Exposure -> M2 + Exposure -> M3 + Exposure -> M4 + M1 -> Outcome + M2 -> Outcome + M3 -> Outcome + M4 -> Outcome + Exposure -> Outcome [style=dashed] +} +") +``` + +In this chapter, we demonstrate multimodal mediation analysis using microbiome +data from the iHMP IBDMDB study from the R/Bioconductor package +curatedMetagenomicData. [@Lloyd-Price2019]. Unlike conventional mediation +analysis, which typically focuses on a single mediator or a small set of +mediators, we performed multimodal mediation analysis with R/Bioconductor +package \[multimedia\] [@Jiang2025] here, which can handle many potential +mediators across multiple data modalities. For example, species-level taxonomic +abundances, functional pathways, and metabolomic profiles can all serve as +potential mediators. This method allows us to estimate and rank each mediator’s +contribution while accounting for correlations between mediators in a +high-dimensional setting. + +In this example, time point (baseline vs. post-treatment) is used as the +treatment, and microbiome-based dysbiosis scores are used as the outcome. We +will first focus on the relative abundances of gut microbial species as +mediators, followed by functional pathways. + +Generally, we will proceed through the following key steps: + +1. Prepare the microbiome taxonomic data from the iHMP IBDMDB project. + +2. Calculate dysbiosis scores as the outcome. + +3. Define the mediation analysis data structure. + +4. Fit the multimodal mediation model using the R package multimedia. + +5. Interpret both the overall indirect effects and the mediator-specific +indirect effects. + +6. Visualize the results using forest plots, histograms, and rankings to +prioritize findings. + +7. Repeat the process for the iHMP microbial pathways. + +### Performing multimodal mediation analysis for iHMP species relative abundance {#sec-relative-abundance} + +We begin by loading the relative abundance data and calculating dysbiosis scores +for each sample based on Bray–Curtis dissimilarity from a healthy reference set. +This represents deviation from a healthy microbiome, where higher dysbiosis +scores indicate greater microbial imbalance. We then subset to IBD subjects with +samples collected at both baseline (visit 1) and post-treatment (visit 21). + +```{r} +################## +# Load libraries # +################## + +library(curatedMetagenomicData) +library(SummarizedExperiment) +library(dplyr) +library(vegan) +library(tidyverse) +library(multimedia) +library(stringr) +library(ggplot2) + + +################## +# Load iHMP data # +################## + +# Load data +se_relative <- curatedMetagenomicData( + "HMP_2019_ibdmdb.relative_abundance", + dryrun = FALSE + )[[1]] + +# Normalize to proportions (from %) +assay(se_relative) <- assay(se_relative) / 100 + +# Change the rownames names for variables of interest +rownames(se_relative) <- sub('.*s__', '', rownames(se_relative)) + + +######################## +# Reference nonIBD set # +######################## + +# Reference set: healthy +ref_set <- colData(se_relative)$disease == "healthy" + + +######################################## +# Calculate Bray-Curtis dissimilarity # +######################################## + +# Bray-Curtis dissimilarity +dist <- as.matrix(vegdist(t(assay(se_relative)), method = "bray", na.rm = TRUE)) + +# Assign SampleID for matching +colData(se_relative)$SampleID <- colnames(se_relative) + + +################################# +# Calculate the dysbiosis score # +################################# + +# Calculate dysbiosis score for each sample +colData(se_relative)$dysbiosis <- sapply(seq_along(ref_set), function(i) { + median(dist[i, ref_set & + (colData(se_relative)$SampleID != colData(se_relative)$SampleID[i])], + na.rm = TRUE) +}) + + +############################################################# +# Subset to IBD only and only one post-baseline time point # +############################################################# + +# Keep IBD only +se_relative <- se_relative[, colData(se_relative)$disease != "healthy"] + +# Visit 1 (baseline) or 25 (post) +se_relative <- se_relative[, colData(se_relative)$visit_number %in% c(1, 21)] + +# Keep subjects with both visits +keep_subjects <- names(which(table(colData(se_relative)$subject_id) == 2)) +se_relative <- se_relative[, colData(se_relative)$subject_id %in% keep_subjects] + +# Define time point label +colData(se_relative)$Time_point <- ifelse( + colData(se_relative)$visit_number == 1, "T0", "T1" +) + +``` + +Next, we will define the mediation anlaysis components. By using the processed +data, we define: + +- Treatment: time point (baseline vs. post-baseline) +- Outcome: dysbiosis score +- Mediators: scaled species abundances + +These are bundled into a mediation data object (exper), which is passed into the +multimedia function to estimate path models for the treatment-to-mediator +(α path), mediator-to-outcome (β path), and the combined indirect effect (α×β +path). + +```{r define-mediation-data-frame} + +############################# +# Define mediation data set # +############################# + +# Treatment: Pre and Post (Baseline vs. Post-baseline) +treatment <- factor(colData(se_relative)$Time_point, levels = c("T0", "T1")) + +# Outcome: Dysbiosis Score +dysbiosis <- as.numeric(colData(se_relative)$dysbiosis) + +# Mediators (scaled, directly from SE): Species +M <- t(assay(se_relative)) +M <- scale(M) +M_clean <- M[, colSums(is.na(M)) == 0] + +# clean mediator names +raw <- colnames(M_clean) +clean1 <- str_remove_all(raw, "\\[|\\]") +clean2 <- str_replace_all(clean1, "[: \\.,]", "_") + +# Save the final safe names +safe_names <- make.names(clean2, unique = TRUE) + +# Apply names +colnames(M_clean) <- safe_names + +# Convert to tibble +mediators <- as_tibble(M_clean) + +# Create the mediation data frame +df <- data.frame( + treatment = treatment, + dysbiosis = dysbiosis +) %>% + bind_cols(mediators) + +``` + +We now fit the multimodal mediation model and inspect the overall indirect and +direct effects to understand whether microbiome composition mediates the +treatment's effect on dysbiosis. + +```{r overall-mediation-analysis} + +############################# +# Run mediation (non-delta) # +############################# + +# Create the Mediation Data object +exper <- mediation_data( + df, + outcome = "dysbiosis", + treatment = "treatment", + mediators = colnames(mediators) +) + +# Fit and summarize +mdl <- multimedia(exper) +res <- estimate(mdl, exper) + +# Summarize overall indirect/direct +summary(res) +print(indirect_overall(res, exper)) +print(direct_effect(res, exper)) + +``` + +```{r specific-mediation} + +############################# +# Mediator-specific effects +############################# + +# Extract the treatment-to-mediator path coefficients +treatment_to_mediator_coef <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) +names(treatment_to_mediator_coef) <- sub("\\.treatmentT1$", "", names(treatment_to_mediator_coef)) + +# Extract the mediator-to-outcome path coefficients (adjusted for all mediators) +mediator_to_outcome_coef <- coef(res@outcome@estimates$dysbiosis)[names(treatment_to_mediator_coef)] + +# Calculate the indirect effect (product of path coefficients) +indirect_effects <- treatment_to_mediator_coef * mediator_to_outcome_coef + +# Summary table +mediator_effects <- data.frame( + mediator = names(treatment_to_mediator_coef), + treatment_to_mediator_coef = treatment_to_mediator_coef, + mediator_to_outcome_coef = mediator_to_outcome_coef, + indirect_effect = indirect_effects +) + +# Rank +top_mediators <- mediator_effects %>% + arrange(desc(abs(indirect_effect))) %>% + head(10) + +print(top_mediators) + +``` + +To quantify uncertainty around the indirect effects, we use non-parametric +bootstrap resampling to compute **95% confidence intervals** for both the +overall and mediator-specific indirect effects. + +```{r visualization-for-overall-indirect-effects} + +##################################### +# Bootstrap overall indirect effect # +##################################### + +set.seed(12345) +boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) +summary(boot_overall$indirect) +quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) + +# Histogram +ggplot(boot_overall$indirect) + + geom_histogram(aes(indirect_effect), bins = 20, fill = "#69b3a2", color = "black") + + theme_classic() + + labs( + x = "Overall Indirect Effect", + y = "Frequency", + title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" + ) + +``` + +Finally, we visualize the **mediator-specific indirect effects** using a forest +plot, displaying point estimates and 95% CIs to highlight the most influential +mediators. + +```{r visualization-for-mediation-specific-indirect-effects} + +############################################### +# Bootstrap mediator-specific indirect effect # +############################################### + +# Define an mediation-specific indirect effect function manually +indirect_each <- function(mdl, exper) { + res <- estimate(mdl, exper) + alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) + names(alpha) <- sub("\\.treatmentT1$", "", names(alpha)) + beta <- coef(res@outcome@estimates$dysbiosis)[names(alpha)] + indirect <- alpha * beta + names(indirect) <- names(alpha) + indirect <- indirect[!is.na(indirect)] + return(indirect) +} + +set.seed(12345) +boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) +summary(boot_each$indirect) + +# Summarize the bootstrap results +lower_upper <- apply(boot_each$indirect, 2, function(x) quantile(x, c(0.025, 0.975), na.rm=TRUE)) +means <- apply(boot_each$indirect, 2, mean, na.rm=TRUE) + +# Convert to data frame +summary_df <- data.frame( + mediator = colnames(boot_each$indirect), + estimate = means, + lower = lower_upper[1,], + upper = lower_upper[2,] +) +summary_df <- summary_df[-1, ] + +# Mark significance if CI does not cross zero +summary_df$significant <- with(summary_df, lower > 0 | upper < 0) + +# Remove the species with 0 confidence interval length +summary_df <- summary_df %>% + filter((upper - lower) != 0) + +# Add rankings +summary_df$rank <- ifelse(summary_df$significant, 1, 2) # 1 for significant, 2 for not +summary_df <- summary_df[order(summary_df$rank, -abs(summary_df$estimate)), ] + +# Forest plot +ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=significant)) + + geom_point() + + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0.2) + + geom_vline(xintercept=0, linetype="dashed", color="grey") + + theme_classic() + + labs( + x = "Observed Indirect Effect with Bootstrap CI", + y = "Mediator", + title = "Forest Plot of Mediation-specific Indirect Effects for Species" + ) + + scale_color_manual(values=c("black","red")) + + theme(legend.position="bottom") + +``` + +Based on the forest plot results, we can see that the agathobaculum +butyriciproducens and escherichia coli pathways has significant indirect +mediation effects. + +Next we repeat the mediation analysis with the iHMP pathway abundances using the +same setup as before, using the time point as the treatment, dysbiosis score as +the outcome, but functional pathway abundance as the mediation variables instead +of the relative species abundance. + +### Performing multimodal mediation analysis for functional pathway abundance {#sec-pathway-abundance} + +Let's extract the pathway abundance data from the iHMP data first, as well as +other data pre-processing. + +```{r load-pkg-data} + +se_pathway <- curatedMetagenomicData( + "HMP_2019_ibdmdb.pathway_abundance", + dryrun = FALSE + )[[1]] + +rows_to_keep <- !grepl("\\|", rownames(se_pathway)) +rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE +se_pathway <- se_pathway[rows_to_keep, ] +assay(se_pathway) <- assay(se_pathway) / 100 +ref_set <- colData(se_pathway)$disease == "healthy" + +dist <- as.matrix(vegdist(t(assay(se_pathway)), method = "bray", na.rm = TRUE)) +colData(se_pathway)$SampleID <- colnames(se_pathway) +colData(se_pathway)$dysbiosis <- sapply(seq_along(ref_set), function(i) { + median(dist[i, ref_set & + (colData(se_pathway)$SampleID != colData(se_pathway)$SampleID[i])], + na.rm = TRUE) +}) + +se_pathway <- se_pathway[, colData(se_pathway)$disease != "healthy"] +se_pathway <- se_pathway[, colData(se_pathway)$visit_number %in% c(1, 27)] +keep_subjects <- names(which(table(colData(se_pathway)$subject_id) == 2)) +se_pathway <- se_pathway[, colData(se_pathway)$subject_id %in% keep_subjects] + +colData(se_pathway)$Time_point <- ifelse( + colData(se_pathway)$visit_number == 1, "T0", "T1" +) + +``` + +We will define the mediation analysis data set then. + +```{r define-mediation-data-frame} + +treatment <- factor(colData(se_pathway)$Time_point, levels = c("T0", "T1")) +dysbiosis <- as.numeric(colData(se_pathway)$dysbiosis) + +M <- t(assay(se_pathway)) +M <- scale(M) +M_clean <- M[, colSums(is.na(M)) == 0] + +raw <- colnames(M_clean) +clean1 <- str_remove_all(raw, "\\[|\\]") +clean2 <- str_replace_all(clean1, "[: \\.,]", "_") +safe_names <- make.names(clean2, unique = TRUE) +colnames(M_clean) <- safe_names + +mediators <- as_tibble(M_clean) + +df <- data.frame( + treatment = treatment, + dysbiosis = dysbiosis +) %>% + bind_cols(mediators) + +``` + +Next, we continue to fit the multimodal mediation model, and extract the overall +and mediation-specific indirect effects. + +```{r mediation-analysis} + +exper <- mediation_data( + df, + outcome = "dysbiosis", + treatment = "treatment", + mediators = colnames(mediators) +) +mdl <- multimedia(exper) +res <- estimate(mdl, exper) +summary(res) +print(indirect_overall(res, exper)) +print(direct_effect(res, exper)) + + +treatment_to_mediator_coef <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) +names(treatment_to_mediator_coef) <- sub("\\.treatmentT1$", "", names(treatment_to_mediator_coef)) +mediator_to_outcome_coef <- coef(res@outcome@estimates$dysbiosis)[names(treatment_to_mediator_coef)] +indirect_effects <- treatment_to_mediator_coef * mediator_to_outcome_coef + +mediator_effects <- data.frame( + mediator = names(treatment_to_mediator_coef), + treatment_to_mediator_coef = treatment_to_mediator_coef, + mediator_to_outcome_coef = mediator_to_outcome_coef, + indirect_effect = indirect_effects +) + +top_mediators <- mediator_effects %>% + arrange(desc(abs(indirect_effect))) %>% + head(10) +print(top_mediators) + +``` + +The non-parametric bootstrap resampling is also used to compute **95% CIs** for +both the overall and mediator-specific pathway indirect effects. + +```{r visualization-for-overall-indirect-effects} + +set.seed(1234) +boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) +summary(boot_overall$indirect) +quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) + +ggplot(boot_overall$indirect) + + geom_histogram(aes(indirect_effect), bins = 20, fill = "#69b3a2", color = "black") + + theme_classic() + + labs( + x = "Overall Indirect Effect", + y = "Frequency", + title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" + ) + +``` + +We visualize the **mediator-specific indirect effects** using a forest plot for +pathway aubundance as well. with the point estimates and 95% CIs. + +```{r visualization-for-mediation-specific-indirect-effects} + +indirect_each <- function(mdl, exper) { + res <- estimate(mdl, exper) + alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) + names(alpha) <- sub("\\.treatmentT1$", "", names(alpha)) + beta <- coef(res@outcome@estimates$dysbiosis)[names(alpha)] + indirect <- alpha * beta + names(indirect) <- names(alpha) + indirect <- indirect[!is.na(indirect)] + return(indirect) +} + +set.seed(1234) +boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) + +lower_upper <- apply(boot_each$indirect, 2, function(x) quantile(x, c(0.025, 0.975), na.rm=TRUE)) +means <- apply(boot_each$indirect, 2, mean, na.rm=TRUE) +summary_df <- data.frame( + mediator = colnames(boot_each$indirect), + estimate = means, + lower = lower_upper[1,], + upper = lower_upper[2,] +) + +summary_df <- summary_df[-1, ] +summary_df$significant <- with(summary_df, lower > 0 | upper < 0) +summary_df <- summary_df[summary_df$upper != summary_df$lower, ] + +# Clean the pathway names +pwy <- rownames(summary_df) +pwy <- gsub("PWY0\\.", "PWY0-", pwy) # PWY0.xxx → PWY0-xxx +pwy <- gsub("PWY\\.", "PWY-", pwy) # PWY.xxx → PWY-xxx +pwy <- gsub("\\.", " ", pwy) # leftover dots → spaces +pwy <- gsub("__", ": ", pwy) # double underscores → colon +pwy <- gsub("_\\.", " ", pwy) # underscore then dot → space +pwy <- gsub("_", " ", pwy) # remaining underscores → space +pwy <- gsub("\\.\\.", " ", pwy) # double dots → space +pwy <- gsub("\\.$", "", pwy) # remove ending period +# pwy <- trimws(pwy) # final cleanup +rownames(summary_df) <- pwy + +ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=significant)) + + geom_point() + + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0.2) + + geom_vline(xintercept=0, linetype="dashed", color="grey") + + theme_classic() + + labs( + x = "Observed Indirect Effect with Bootstrap CI", + y = "Mediator", + title = "Forest Plot of Mediation-specific Indirect Effects for Pathway" + ) + + scale_color_manual(values=c("black","red")) + + theme(legend.position="bottom") + +``` + +Based on the forest plot results, we can see that the tRNA charging pathway has +significant indirect mediation effects on the dysbiosis scores. From af9fafef80726ba50c8b6de03cc2d5dd44973f7d Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Thu, 7 Aug 2025 04:04:51 +0800 Subject: [PATCH 02/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 223 +++++++++++++++++++++++--------------- 1 file changed, 135 insertions(+), 88 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 439b2666..c17292e5 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -7,7 +7,10 @@ author: Yihan Liu, Himel Mallick # **Multimodal Mediation Analysis** {#sec-MSEA} -```{r setup, echo=FALSE, results="asis"} +```{r} +#| label: setup +#| echo: false +#| results: asis remove(list = ls()) invisible(gc()) library(rebook) @@ -28,7 +31,7 @@ layers, providing mechanistic insights into **how exposures impact outcomes** **via molecular pathways**. ```{r} -#| label: fig-multimodal-mediation +#| label: fig_multimodal_mediation #| fig-cap: Directed acyclic graph illustrating multimodal mediation where an exposure affects multiple mediators across modalities (microbial species and metabolites), which in turn affect the outcome. A direct path from exposure to outcome is also included (dashed). While gut microbial species biologically influence metabolite levels, for tractability in high-dimensional mediation analysis we assume species and metabolites act as conditionally independent parallel mediators given exposure, acknowledging this as an approximation #| echo: false @@ -68,7 +71,7 @@ analysis, which typically focuses on a single mediator or a small set of mediators, we performed multimodal mediation analysis with R/Bioconductor package \[multimedia\] [@Jiang2025] here, which can handle many potential mediators across multiple data modalities. For example, species-level taxonomic -abundances, functional pathways, and metabolomic profiles can all serve as +abundances, pathways abundances, and metabolomic profiles can all serve as potential mediators. This method allows us to estimate and rank each mediator’s contribution while accounting for correlations between mediators in a high-dimensional setting. @@ -76,7 +79,7 @@ high-dimensional setting. In this example, time point (baseline vs. post-treatment) is used as the treatment, and microbiome-based dysbiosis scores are used as the outcome. We will first focus on the relative abundances of gut microbial species as -mediators, followed by functional pathways. +mediators, followed by pathways abundances. Generally, we will proceed through the following key steps: @@ -105,6 +108,9 @@ scores indicate greater microbial imbalance. We then subset to IBD subjects with samples collected at both baseline (visit 1) and post-treatment (visit 21). ```{r} +#| label: iHMP_data_preprocessing +#| message: false + ################## # Load libraries # ################## @@ -117,78 +123,94 @@ library(tidyverse) library(multimedia) library(stringr) library(ggplot2) - +library(mia) ################## # Load iHMP data # ################## # Load data -se_relative <- curatedMetagenomicData( +tse_relative <- curatedMetagenomicData( "HMP_2019_ibdmdb.relative_abundance", dryrun = FALSE )[[1]] -# Normalize to proportions (from %) -assay(se_relative) <- assay(se_relative) / 100 +# Assign SampleID for matching +colData(tse_relative)$SampleID <- colnames(tse_relative) -# Change the rownames names for variables of interest -rownames(se_relative) <- sub('.*s__', '', rownames(se_relative)) +# Convert relative_abundance assay to relabundance (which is in [0,1] interval) +tse_relative <- transformAssay(tse_relative, assay.type="relative_abundance", method="relabundance") + +# Optionally, remove the original assay to avoid confusion between the two relative abundance versions +assay(tse_relative, "relative_abundance") <- NULL +# Remove samples with NA in the relabundance assay +tse_relative <- tse_relative[, colSums(is.na(assay(tse_relative))) == 0] + +# Change the rownames names for variables of interest +rownames(tse_relative) <- sub('.*s__', '', rownames(tse_relative)) +rownames(tse_relative) <- str_remove_all(rownames(tse_relative), "\\[|\\]") +rownames(tse_relative) <- str_replace_all(rownames(tse_relative), "[: \\.,]", "_") +safe_names <- make.names(rownames(tse_relative), unique = TRUE) +rownames(tse_relative) <- safe_names ######################## # Reference nonIBD set # ######################## # Reference set: healthy -ref_set <- colData(se_relative)$disease == "healthy" - +tse_relative$disease_binary <- tse_relative$disease == "healthy" ######################################## # Calculate Bray-Curtis dissimilarity # ######################################## # Bray-Curtis dissimilarity -dist <- as.matrix(vegdist(t(assay(se_relative)), method = "bray", na.rm = TRUE)) - -# Assign SampleID for matching -colData(se_relative)$SampleID <- colnames(se_relative) - +diss <- as.matrix(getDissimilarity(tse_relative, method = "bray", na.rm=TRUE, assay.type = "relabundance")) ################################# # Calculate the dysbiosis score # ################################# # Calculate dysbiosis score for each sample +# For each sample i, we compute the median distance between sample i and all reference samples in `ref_set` +sample_ids <- colData(se_relative)$SampleID + colData(se_relative)$dysbiosis <- sapply(seq_along(ref_set), function(i) { - median(dist[i, ref_set & - (colData(se_relative)$SampleID != colData(se_relative)$SampleID[i])], - na.rm = TRUE) -}) + # Logical vector indicating all other reference samples (excluding i) + ref_others <- ref_set & (sample_ids != sample_ids[i]) + + # Compute median distance between sample i and these reference samples + median(dist[i, ref_others], na.rm = TRUE) + +}) ############################################################# # Subset to IBD only and only one post-baseline time point # ############################################################# # Keep IBD only -se_relative <- se_relative[, colData(se_relative)$disease != "healthy"] +tse_relative <- tse_relative[, colData(tse_relative)$disease != "healthy"] # Visit 1 (baseline) or 25 (post) -se_relative <- se_relative[, colData(se_relative)$visit_number %in% c(1, 21)] +tse_relative <- tse_relative[, colData(tse_relative)$visit_number %in% c(1, 21)] # Keep subjects with both visits -keep_subjects <- names(which(table(colData(se_relative)$subject_id) == 2)) -se_relative <- se_relative[, colData(se_relative)$subject_id %in% keep_subjects] +keep_subjects <- names(which(table(colData(tse_relative)$subject_id) == 2)) +tse_relative <- tse_relative[, colData(tse_relative)$subject_id %in% keep_subjects] # Define time point label -colData(se_relative)$Time_point <- ifelse( - colData(se_relative)$visit_number == 1, "T0", "T1" +colData(tse_relative)$Time_point <- ifelse( + colData(tse_relative)$visit_number == 1, "T0", "T1" ) +# Define treatment: Pre and Post (Baseline vs. Post-baseline) +treatment <- factor(colData(tse_relative)$Time_point, levels = c("T0", "T1")) + ``` -Next, we will define the mediation anlaysis components. By using the processed +Next, we will define the mediation analysis components. By using the processed data, we define: - Treatment: time point (baseline vs. post-baseline) @@ -200,41 +222,28 @@ multimedia function to estimate path models for the treatment-to-mediator (α path), mediator-to-outcome (β path), and the combined indirect effect (α×β path). -```{r define-mediation-data-frame} +```{r} +#| label: define_mediation_data_frame +#| message: false ############################# # Define mediation data set # ############################# -# Treatment: Pre and Post (Baseline vs. Post-baseline) -treatment <- factor(colData(se_relative)$Time_point, levels = c("T0", "T1")) - # Outcome: Dysbiosis Score -dysbiosis <- as.numeric(colData(se_relative)$dysbiosis) +colData(tse_relative)$dysbiosis <- as.numeric(colData(tse_relative)$dysbiosis) # Mediators (scaled, directly from SE): Species -M <- t(assay(se_relative)) -M <- scale(M) -M_clean <- M[, colSums(is.na(M)) == 0] - -# clean mediator names -raw <- colnames(M_clean) -clean1 <- str_remove_all(raw, "\\[|\\]") -clean2 <- str_replace_all(clean1, "[: \\.,]", "_") - -# Save the final safe names -safe_names <- make.names(clean2, unique = TRUE) - -# Apply names -colnames(M_clean) <- safe_names +tse_relative <- transformAssay(tse_relative, assay.type="relabundance", method="standardize", name="scaled") +M <- assay(tse_relative, "scaled") # Convert to tibble -mediators <- as_tibble(M_clean) +mediators <- as_tibble(M) # Create the mediation data frame df <- data.frame( treatment = treatment, - dysbiosis = dysbiosis + dysbiosis = tse_relative$dysbiosis ) %>% bind_cols(mediators) @@ -244,7 +253,9 @@ We now fit the multimodal mediation model and inspect the overall indirect and direct effects to understand whether microbiome composition mediates the treatment's effect on dysbiosis. -```{r overall-mediation-analysis} +```{r} +#| label: overall_mediation_analysis +#| message: false ############################# # Run mediation (non-delta) # @@ -269,7 +280,9 @@ print(direct_effect(res, exper)) ``` -```{r specific-mediation} +```{r} +#| label: specific_mediation +#| message: false ############################# # Mediator-specific effects @@ -306,7 +319,9 @@ To quantify uncertainty around the indirect effects, we use non-parametric bootstrap resampling to compute **95% confidence intervals** for both the overall and mediator-specific indirect effects. -```{r visualization-for-overall-indirect-effects} +```{r} +#| label: visualization_for_overall_indirect_effects +#| message: false ##################################### # Bootstrap overall indirect effect # @@ -333,7 +348,9 @@ Finally, we visualize the **mediator-specific indirect effects** using a forest plot, displaying point estimates and 95% CIs to highlight the most influential mediators. -```{r visualization-for-mediation-specific-indirect-effects} +```{r} +#| label: visualization_for_mediation_specific_indirect_effects +#| message: false ############################################### # Bootstrap mediator-specific indirect effect # @@ -395,74 +412,93 @@ ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=signific ``` -Based on the forest plot results, we can see that the agathobaculum -butyriciproducens and escherichia coli pathways has significant indirect +Based on the forest plot results, we can see that the Agathobaculum +butyriciproducens and Escherichia coli pathways has significant indirect mediation effects. Next we repeat the mediation analysis with the iHMP pathway abundances using the same setup as before, using the time point as the treatment, dysbiosis score as -the outcome, but functional pathway abundance as the mediation variables instead +the outcome, but pathways abundances as the mediation variables instead of the relative species abundance. -### Performing multimodal mediation analysis for functional pathway abundance {#sec-pathway-abundance} +### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} Let's extract the pathway abundance data from the iHMP data first, as well as other data pre-processing. -```{r load-pkg-data} +```{r} +#| label: load_pkg_data +#| message: false -se_pathway <- curatedMetagenomicData( +tse_pathway <- curatedMetagenomicData( "HMP_2019_ibdmdb.pathway_abundance", dryrun = FALSE )[[1]] -rows_to_keep <- !grepl("\\|", rownames(se_pathway)) +colData(tse_pathway)$SampleID <- colnames(tse_pathway) +tse_pathway <- tse_pathway[, colSums(is.na(assay(tse_pathway))) == 0] + +# Clean the pathway names +pwy <- rownames(tse_pathway) +pwy <- gsub("PWY0\\.", "PWY0-", pwy) # PWY0.xxx → PWY0-xxx +pwy <- gsub("PWY\\.", "PWY-", pwy) # PWY.xxx → PWY-xxx +pwy <- gsub("\\.", " ", pwy) # leftover dots → spaces +pwy <- gsub("__", ": ", pwy) # double underscores → colon +pwy <- gsub("_\\.", " ", pwy) # underscore then dot → space +pwy <- gsub("_", " ", pwy) # remaining underscores → space +pwy <- gsub("\\.\\.", " ", pwy) # double dots → space +pwy <- gsub("\\.$", "", pwy) # remove ending period +# pwy <- trimws(pwy) # final cleanup +rownames(summary_df) <- pwy + +rows_to_keep <- !grepl("\\|", rownames(tse_pathway)) rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE -se_pathway <- se_pathway[rows_to_keep, ] -assay(se_pathway) <- assay(se_pathway) / 100 -ref_set <- colData(se_pathway)$disease == "healthy" - -dist <- as.matrix(vegdist(t(assay(se_pathway)), method = "bray", na.rm = TRUE)) -colData(se_pathway)$SampleID <- colnames(se_pathway) -colData(se_pathway)$dysbiosis <- sapply(seq_along(ref_set), function(i) { - median(dist[i, ref_set & - (colData(se_pathway)$SampleID != colData(se_pathway)$SampleID[i])], +tse_pathway <- tse_pathway[rows_to_keep, ] +assay(tse_pathway) <- assay(tse_pathway) / 100 +tse_pathway$disease_binary <- colData(tse_pathway)$disease == "healthy" + +diss <- as.matrix(vegdist(t(assay(tse_pathway)), method = "bray", na.rm = TRUE)) +colData(tse_pathway)$dysbiosis <- sapply(seq_along(tse_pathway$disease_binary), function(i) { + median(diss[i, tse_pathway$disease_binary & + (colData(tse_pathway)$SampleID != colData(tse_pathway)$SampleID[i])], na.rm = TRUE) }) -se_pathway <- se_pathway[, colData(se_pathway)$disease != "healthy"] -se_pathway <- se_pathway[, colData(se_pathway)$visit_number %in% c(1, 27)] -keep_subjects <- names(which(table(colData(se_pathway)$subject_id) == 2)) -se_pathway <- se_pathway[, colData(se_pathway)$subject_id %in% keep_subjects] +tse_pathway <- tse_pathway[, colData(tse_pathway)$disease != "healthy"] +tse_pathway <- tse_pathway[, colData(tse_pathway)$visit_number %in% c(1, 27)] +keep_subjects <- names(which(table(colData(tse_pathway)$subject_id) == 2)) +tse_pathway <- tse_pathway[, colData(tse_pathway)$subject_id %in% keep_subjects] -colData(se_pathway)$Time_point <- ifelse( - colData(se_pathway)$visit_number == 1, "T0", "T1" +colData(tse_pathway)$Time_point <- ifelse( + colData(tse_pathway)$visit_number == 1, "T0", "T1" ) +treatment <- factor(colData(tse_pathway)$Time_point, levels = c("T0", "T1")) + ``` We will define the mediation analysis data set then. -```{r define-mediation-data-frame} +```{r} +#| label: define_mediation_data_frame +#| message: false -treatment <- factor(colData(se_pathway)$Time_point, levels = c("T0", "T1")) -dysbiosis <- as.numeric(colData(se_pathway)$dysbiosis) +colData(tse_pathway)$dysbiosis <- as.numeric(colData(tse_pathway)$dysbiosis) -M <- t(assay(se_pathway)) -M <- scale(M) -M_clean <- M[, colSums(is.na(M)) == 0] +tse_pathway <- transformAssay(tse_pathway, assay.type="relabundance", method="standardize", name="scaled") +M <- assay(tse_pathway, "scaled") -raw <- colnames(M_clean) +raw <- colnames(M) clean1 <- str_remove_all(raw, "\\[|\\]") clean2 <- str_replace_all(clean1, "[: \\.,]", "_") safe_names <- make.names(clean2, unique = TRUE) -colnames(M_clean) <- safe_names +colnames(M) <- safe_names -mediators <- as_tibble(M_clean) +mediators <- as_tibble(M) df <- data.frame( treatment = treatment, - dysbiosis = dysbiosis + dysbiosis = tse_pathway$dysbiosis ) %>% bind_cols(mediators) @@ -471,7 +507,9 @@ df <- data.frame( Next, we continue to fit the multimodal mediation model, and extract the overall and mediation-specific indirect effects. -```{r mediation-analysis} +```{r} +#| label: mediation_analysis +#| message: false exper <- mediation_data( df, @@ -508,7 +546,9 @@ print(top_mediators) The non-parametric bootstrap resampling is also used to compute **95% CIs** for both the overall and mediator-specific pathway indirect effects. -```{r visualization-for-overall-indirect-effects} +```{r} +#| label: visualization_for_overall_indirect_effects +#| message: false set.seed(1234) boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) @@ -529,7 +569,9 @@ ggplot(boot_overall$indirect) + We visualize the **mediator-specific indirect effects** using a forest plot for pathway aubundance as well. with the point estimates and 95% CIs. -```{r visualization-for-mediation-specific-indirect-effects} +```{r} +#| label: visualization_for_mediation_specific_indirect_effects +#| message: false indirect_each <- function(mdl, exper) { res <- estimate(mdl, exper) @@ -571,6 +613,11 @@ pwy <- gsub("\\.$", "", pwy) # remove ending period # pwy <- trimws(pwy) # final cleanup rownames(summary_df) <- pwy + + + + + ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=significant)) + geom_point() + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0.2) + From fb8c08057562ade8f6aebcbe0f5ada05ff0bde94 Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Sat, 30 Aug 2025 02:36:57 +0800 Subject: [PATCH 03/16] Update multimedia.qmd Just fixed some minor inconsistence in the name of variables --- inst/pages/multimedia.qmd | 27 +++++++-------------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index c17292e5..40c43a98 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -174,15 +174,15 @@ diss <- as.matrix(getDissimilarity(tse_relative, method = "bray", na.rm=TRUE, as # Calculate dysbiosis score for each sample # For each sample i, we compute the median distance between sample i and all reference samples in `ref_set` -sample_ids <- colData(se_relative)$SampleID +sample_ids <- colData(tse_relative)$SampleID -colData(se_relative)$dysbiosis <- sapply(seq_along(ref_set), function(i) { +colData(tse_relative)$dysbiosis <- sapply(seq_along(tse_relative$disease_binary), function(i) { # Logical vector indicating all other reference samples (excluding i) - ref_others <- ref_set & (sample_ids != sample_ids[i]) + ref_others <- tse_relative$disease_binary & (sample_ids != sample_ids[i]) # Compute median distance between sample i and these reference samples - median(dist[i, ref_others], na.rm = TRUE) + median(diss[i, ref_others], na.rm = TRUE) }) @@ -235,7 +235,7 @@ colData(tse_relative)$dysbiosis <- as.numeric(colData(tse_relative)$dysbiosis) # Mediators (scaled, directly from SE): Species tse_relative <- transformAssay(tse_relative, assay.type="relabundance", method="standardize", name="scaled") -M <- assay(tse_relative, "scaled") +M <- t(assay(tse_relative, "scaled")) # Convert to tibble mediators <- as_tibble(M) @@ -438,19 +438,6 @@ tse_pathway <- curatedMetagenomicData( colData(tse_pathway)$SampleID <- colnames(tse_pathway) tse_pathway <- tse_pathway[, colSums(is.na(assay(tse_pathway))) == 0] -# Clean the pathway names -pwy <- rownames(tse_pathway) -pwy <- gsub("PWY0\\.", "PWY0-", pwy) # PWY0.xxx → PWY0-xxx -pwy <- gsub("PWY\\.", "PWY-", pwy) # PWY.xxx → PWY-xxx -pwy <- gsub("\\.", " ", pwy) # leftover dots → spaces -pwy <- gsub("__", ": ", pwy) # double underscores → colon -pwy <- gsub("_\\.", " ", pwy) # underscore then dot → space -pwy <- gsub("_", " ", pwy) # remaining underscores → space -pwy <- gsub("\\.\\.", " ", pwy) # double dots → space -pwy <- gsub("\\.$", "", pwy) # remove ending period -# pwy <- trimws(pwy) # final cleanup -rownames(summary_df) <- pwy - rows_to_keep <- !grepl("\\|", rownames(tse_pathway)) rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE tse_pathway <- tse_pathway[rows_to_keep, ] @@ -485,8 +472,8 @@ We will define the mediation analysis data set then. colData(tse_pathway)$dysbiosis <- as.numeric(colData(tse_pathway)$dysbiosis) -tse_pathway <- transformAssay(tse_pathway, assay.type="relabundance", method="standardize", name="scaled") -M <- assay(tse_pathway, "scaled") +tse_pathway <- transformAssay(tse_pathway, assay.type="pathway_abundance", method="standardize", name="scaled") +M <- t(assay(tse_pathway, "scaled")) raw <- colnames(M) clean1 <- str_remove_all(raw, "\\[|\\]") From 9981a192b46935e2fe6069b7fe8406f97c2563c0 Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Fri, 19 Sep 2025 06:20:46 +0800 Subject: [PATCH 04/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 224 ++++++++++++++------------------------ 1 file changed, 83 insertions(+), 141 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 40c43a98..09cd2dda 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -18,17 +18,7 @@ chapterPreamble() ``` -Building upon the concepts and workflows presented in the previous chapter on -mediation analysis, here we demonstrate how to perform **multimodal mediation** -**analysis** to identify potential mediators across multiple omics layers. -**Multimodal mediation analysis** examines whether and how an exposure (e.g., -treatment, diet) affects an outcome (e.g., disease, phenotype) **through** -**intermediate variables (mediators) across multiple omics layers** -**simultaneously** (e.g., microbiome, metabolomics). It identifies which -features across modalities act as mediators, quantifies indirect (mediated) -effects while adjusting for covariates, and compares mediation strength across -layers, providing mechanistic insights into **how exposures impact outcomes** -**via molecular pathways**. +Building upon the concepts and workflows presented in the previous chapter on mediation analysis, here we demonstrate how to perform **multimodal mediation** **analysis** to identify potential mediators across multiple omics layers. **Multimodal mediation analysis** examines whether and how an exposure (e.g., treatment, diet) affects an outcome (e.g., disease, phenotype) **through** **intermediate variables (mediators) across multiple omics layers** **simultaneously** (e.g., microbiome, metabolomics). It identifies which features across modalities act as mediators, quantifies indirect (mediated) effects while adjusting for covariates, and compares mediation strength across layers, providing mechanistic insights into **how exposures impact outcomes** **via molecular pathways**. ```{r} #| label: fig_multimodal_mediation @@ -64,22 +54,9 @@ digraph multimodal_mediation { ") ``` -In this chapter, we demonstrate multimodal mediation analysis using microbiome -data from the iHMP IBDMDB study from the R/Bioconductor package -curatedMetagenomicData. [@Lloyd-Price2019]. Unlike conventional mediation -analysis, which typically focuses on a single mediator or a small set of -mediators, we performed multimodal mediation analysis with R/Bioconductor -package \[multimedia\] [@Jiang2025] here, which can handle many potential -mediators across multiple data modalities. For example, species-level taxonomic -abundances, pathways abundances, and metabolomic profiles can all serve as -potential mediators. This method allows us to estimate and rank each mediator’s -contribution while accounting for correlations between mediators in a -high-dimensional setting. - -In this example, time point (baseline vs. post-treatment) is used as the -treatment, and microbiome-based dysbiosis scores are used as the outcome. We -will first focus on the relative abundances of gut microbial species as -mediators, followed by pathways abundances. +In this chapter, we demonstrate multimodal mediation analysis using microbiome data from the iHMP IBDMDB study from the R/Bioconductor package curatedMetagenomicData. [@Lloyd-Price2019]. Unlike conventional mediation analysis, which typically focuses on a single mediator or a small set of mediators, we performed multimodal mediation analysis with R/Bioconductor package \[multimedia\] [@Jiang2025] here, which can handle many potential mediators across multiple data modalities. For example, species-level taxonomic abundances, pathways abundances, and metabolomic profiles can all serve as potential mediators. This method allows us to estimate and rank each mediator’s contribution while accounting for correlations between mediators in a high-dimensional setting. + +In this example, time point (baseline vs. post-treatment) is used as the treatment, and microbiome-based dysbiosis scores are used as the outcome. We will first focus on the relative abundances of gut microbial species as mediators, followed by pathways abundances. Generally, we will proceed through the following key steps: @@ -91,21 +68,15 @@ Generally, we will proceed through the following key steps: 4. Fit the multimodal mediation model using the R package multimedia. -5. Interpret both the overall indirect effects and the mediator-specific -indirect effects. +5. Interpret both the overall indirect effects and the mediator-specific indirect effects. -6. Visualize the results using forest plots, histograms, and rankings to -prioritize findings. +6. Visualize the results using forest plots, histograms, and rankings to prioritize findings. 7. Repeat the process for the iHMP microbial pathways. ### Performing multimodal mediation analysis for iHMP species relative abundance {#sec-relative-abundance} -We begin by loading the relative abundance data and calculating dysbiosis scores -for each sample based on Bray–Curtis dissimilarity from a healthy reference set. -This represents deviation from a healthy microbiome, where higher dysbiosis -scores indicate greater microbial imbalance. We then subset to IBD subjects with -samples collected at both baseline (visit 1) and post-treatment (visit 21). +We begin by loading the relative abundance data and calculating dysbiosis scores for each sample based on Bray–Curtis dissimilarity from a healthy reference set. This represents deviation from a healthy microbiome, where higher dysbiosis scores indicate greater microbial imbalance. We then subset to IBD subjects with samples collected at both baseline (visit 1) and post-treatment (visit 21). ```{r} #| label: iHMP_data_preprocessing @@ -114,7 +85,8 @@ samples collected at both baseline (visit 1) and post-treatment (visit 21). ################## # Load libraries # ################## - +remove(list = ls()) +gc() library(curatedMetagenomicData) library(SummarizedExperiment) library(dplyr) @@ -206,21 +178,17 @@ colData(tse_relative)$Time_point <- ifelse( ) # Define treatment: Pre and Post (Baseline vs. Post-baseline) -treatment <- factor(colData(tse_relative)$Time_point, levels = c("T0", "T1")) +colData(tse_relative)$treatment <- factor(colData(tse_relative)$Time_point, levels = c("T0", "T1")) ``` -Next, we will define the mediation analysis components. By using the processed -data, we define: +Next, we will define the mediation analysis components. By using the processed data, we define: - Treatment: time point (baseline vs. post-baseline) - Outcome: dysbiosis score - Mediators: scaled species abundances -These are bundled into a mediation data object (exper), which is passed into the -multimedia function to estimate path models for the treatment-to-mediator -(α path), mediator-to-outcome (β path), and the combined indirect effect (α×β -path). +These are bundled into a mediation data object (exper), which is passed into the multimedia function to estimate path models for the treatment-to-mediator (α path), mediator-to-outcome (β path), and the combined indirect effect (α×β path). ```{r} #| label: define_mediation_data_frame @@ -235,23 +203,11 @@ colData(tse_relative)$dysbiosis <- as.numeric(colData(tse_relative)$dysbiosis) # Mediators (scaled, directly from SE): Species tse_relative <- transformAssay(tse_relative, assay.type="relabundance", method="standardize", name="scaled") -M <- t(assay(tse_relative, "scaled")) - -# Convert to tibble -mediators <- as_tibble(M) - -# Create the mediation data frame -df <- data.frame( - treatment = treatment, - dysbiosis = tse_relative$dysbiosis -) %>% - bind_cols(mediators) +assays(tse_relative) <- SimpleList(mediators = assay(tse_relative, "scaled")) ``` -We now fit the multimodal mediation model and inspect the overall indirect and -direct effects to understand whether microbiome composition mediates the -treatment's effect on dysbiosis. +We now fit the multimodal mediation model and inspect the overall indirect and direct effects to understand whether microbiome composition mediates the treatment's effect on dysbiosis. ```{r} #| label: overall_mediation_analysis @@ -261,12 +217,24 @@ treatment's effect on dysbiosis. # Run mediation (non-delta) # ############################# +# Convert the TreeSummarizedExperiment to be SummarizedExperiment +se_relative <- tryCatch(as(tse_relative, "SummarizedExperiment"), error = function(e) NULL) +if (is.null(rownames(se_relative)) || any(rownames(se_relative) == "")) { + rn <- rownames(tse_relative) + stopifnot(length(rn) == nrow(se_relative)) + rownames(se_relative) <- rn + rownames(rowData(se_relative)) <- rn # keep rowData aligned too +} + +# Use indices for mediators +medi_idx <- seq_len(nrow(se_relative)) + # Create the Mediation Data object -exper <- mediation_data( - df, - outcome = "dysbiosis", - treatment = "treatment", - mediators = colnames(mediators) +exper <- multimedia::mediation_data( + se_relative, + outcomes = "dysbiosis", + treatments = "treatment", + mediators = medi_idx ) # Fit and summarize @@ -280,6 +248,7 @@ print(direct_effect(res, exper)) ``` + ```{r} #| label: specific_mediation #| message: false @@ -288,26 +257,27 @@ print(direct_effect(res, exper)) # Mediator-specific effects ############################# -# Extract the treatment-to-mediator path coefficients -treatment_to_mediator_coef <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) -names(treatment_to_mediator_coef) <- sub("\\.treatmentT1$", "", names(treatment_to_mediator_coef)) +# Helper function to extract the effects +extract_effects <- function(res, treat_term = "treatmentT1", outcome = "dysbiosis") { + # α: treatment -> each mediator (one model per mediator) + alpha <- sapply(res@mediation@estimates, function(m) unname(coef(m)[treat_term])) -# Extract the mediator-to-outcome path coefficients (adjusted for all mediators) -mediator_to_outcome_coef <- coef(res@outcome@estimates$dysbiosis)[names(treatment_to_mediator_coef)] + # β: mediator -> outcome (single model with all mediators) + beta <- unname(coef(res@outcome@estimates[[outcome]])[names(alpha)]) -# Calculate the indirect effect (product of path coefficients) -indirect_effects <- treatment_to_mediator_coef * mediator_to_outcome_coef + tibble::tibble( + mediator = names(alpha), + alpha = as.numeric(alpha), + beta = as.numeric(beta), + indirect_effect = alpha * beta + ) +} -# Summary table -mediator_effects <- data.frame( - mediator = names(treatment_to_mediator_coef), - treatment_to_mediator_coef = treatment_to_mediator_coef, - mediator_to_outcome_coef = mediator_to_outcome_coef, - indirect_effect = indirect_effects -) +# Extract the treatment-to-mediator path coefficients +effects_by_mediator <- extract_effects(res, treat_term = "treatmentT1", outcome = "dysbiosis") # Rank -top_mediators <- mediator_effects %>% +top_mediators <- effects_by_mediator %>% arrange(desc(abs(indirect_effect))) %>% head(10) @@ -315,9 +285,7 @@ print(top_mediators) ``` -To quantify uncertainty around the indirect effects, we use non-parametric -bootstrap resampling to compute **95% confidence intervals** for both the -overall and mediator-specific indirect effects. +To quantify uncertainty around the indirect effects, we use non-parametric bootstrap resampling to compute **95% confidence intervals** for both the overall and mediator-specific indirect effects. ```{r} #| label: visualization_for_overall_indirect_effects @@ -344,9 +312,7 @@ ggplot(boot_overall$indirect) + ``` -Finally, we visualize the **mediator-specific indirect effects** using a forest -plot, displaying point estimates and 95% CIs to highlight the most influential -mediators. +Finally, we visualize the **mediator-specific indirect effects** using a forest plot, displaying point estimates and 95% CIs to highlight the most influential mediators. ```{r} #| label: visualization_for_mediation_specific_indirect_effects @@ -412,19 +378,13 @@ ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=signific ``` -Based on the forest plot results, we can see that the Agathobaculum -butyriciproducens and Escherichia coli pathways has significant indirect -mediation effects. +Based on the forest plot results, we can see that the Agathobaculum butyriciproducens and Escherichia coli pathways has significant indirect mediation effects. -Next we repeat the mediation analysis with the iHMP pathway abundances using the -same setup as before, using the time point as the treatment, dysbiosis score as -the outcome, but pathways abundances as the mediation variables instead -of the relative species abundance. +Next we repeat the mediation analysis with the iHMP pathway abundances using the same setup as before, using the time point as the treatment, dysbiosis score as the outcome, but pathways abundances as the mediation variables instead of the relative species abundance. ### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} -Let's extract the pathway abundance data from the iHMP data first, as well as -other data pre-processing. +Let's extract the pathway abundance data from the iHMP data first, as well as other data pre-processing. ```{r} #| label: load_pkg_data @@ -460,7 +420,7 @@ colData(tse_pathway)$Time_point <- ifelse( colData(tse_pathway)$visit_number == 1, "T0", "T1" ) -treatment <- factor(colData(tse_pathway)$Time_point, levels = c("T0", "T1")) +colData(tse_pathway)$treatment <- factor(colData(tse_pathway)$Time_point, levels = c("T0", "T1")) ``` @@ -471,67 +431,56 @@ We will define the mediation analysis data set then. #| message: false colData(tse_pathway)$dysbiosis <- as.numeric(colData(tse_pathway)$dysbiosis) - tse_pathway <- transformAssay(tse_pathway, assay.type="pathway_abundance", method="standardize", name="scaled") -M <- t(assay(tse_pathway, "scaled")) - -raw <- colnames(M) -clean1 <- str_remove_all(raw, "\\[|\\]") -clean2 <- str_replace_all(clean1, "[: \\.,]", "_") -safe_names <- make.names(clean2, unique = TRUE) -colnames(M) <- safe_names - -mediators <- as_tibble(M) - -df <- data.frame( - treatment = treatment, - dysbiosis = tse_pathway$dysbiosis -) %>% - bind_cols(mediators) +assays(tse_pathway) <- SimpleList(mediators = assay(tse_pathway, "scaled")) ``` -Next, we continue to fit the multimodal mediation model, and extract the overall -and mediation-specific indirect effects. +Next, we continue to fit the multimodal mediation model, and extract the overall and mediation-specific indirect effects. ```{r} #| label: mediation_analysis #| message: false -exper <- mediation_data( - df, - outcome = "dysbiosis", - treatment = "treatment", - mediators = colnames(mediators) +se_pathway <- tryCatch(as(tse_pathway, "SummarizedExperiment"), error = function(e) NULL) +if (is.null(rownames(se_pathway)) || any(rownames(se_pathway) == "")) { + rn <- rownames(tse_pathway) + stopifnot(length(rn) == nrow(se_pathway)) + rownames(se_pathway) <- rn + rownames(rowData(se_pathway)) <- rn # keep rowData aligned too +} + +raw <- rownames(se_pathway) +clean1 <- str_remove_all(raw, "\\[|\\]") +clean2 <- str_replace_all(clean1, "[: \\.,]", "_") +safe_names <- make.names(clean2, unique = TRUE) +rownames(se_pathway) <- safe_names + +medi_idx <- seq_len(nrow(se_pathway)) + +exper <- multimedia::mediation_data( + se_pathway, + outcomes = "dysbiosis", + treatments = "treatment", + mediators = medi_idx ) + mdl <- multimedia(exper) res <- estimate(mdl, exper) summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) +effects_by_mediator <- extract_effects(res, treat_term = "treatmentT1", outcome = "dysbiosis") -treatment_to_mediator_coef <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) -names(treatment_to_mediator_coef) <- sub("\\.treatmentT1$", "", names(treatment_to_mediator_coef)) -mediator_to_outcome_coef <- coef(res@outcome@estimates$dysbiosis)[names(treatment_to_mediator_coef)] -indirect_effects <- treatment_to_mediator_coef * mediator_to_outcome_coef - -mediator_effects <- data.frame( - mediator = names(treatment_to_mediator_coef), - treatment_to_mediator_coef = treatment_to_mediator_coef, - mediator_to_outcome_coef = mediator_to_outcome_coef, - indirect_effect = indirect_effects -) - -top_mediators <- mediator_effects %>% +top_mediators <- effects_by_mediator %>% arrange(desc(abs(indirect_effect))) %>% head(10) print(top_mediators) ``` -The non-parametric bootstrap resampling is also used to compute **95% CIs** for -both the overall and mediator-specific pathway indirect effects. +The non-parametric bootstrap resampling is also used to compute **95% CIs** for both the overall and mediator-specific pathway indirect effects. ```{r} #| label: visualization_for_overall_indirect_effects @@ -553,8 +502,7 @@ ggplot(boot_overall$indirect) + ``` -We visualize the **mediator-specific indirect effects** using a forest plot for -pathway aubundance as well. with the point estimates and 95% CIs. +We visualize the **mediator-specific indirect effects** using a forest plot for pathway aubundance as well. with the point estimates and 95% CIs. ```{r} #| label: visualization_for_mediation_specific_indirect_effects @@ -600,11 +548,6 @@ pwy <- gsub("\\.$", "", pwy) # remove ending period # pwy <- trimws(pwy) # final cleanup rownames(summary_df) <- pwy - - - - - ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=significant)) + geom_point() + geom_errorbarh(aes(xmin=lower, xmax=upper), height=0.2) + @@ -620,5 +563,4 @@ ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=signific ``` -Based on the forest plot results, we can see that the tRNA charging pathway has -significant indirect mediation effects on the dysbiosis scores. +Based on the forest plot results, we can see that the tRNA charging pathway has significant indirect mediation effects on the dysbiosis scores. From d2be9147566822bbe24e998dd87135fed3ae6dac Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Fri, 21 Nov 2025 07:37:27 +0800 Subject: [PATCH 05/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 211 +++++++++++++++++++++----------------- 1 file changed, 117 insertions(+), 94 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 09cd2dda..18ef63f1 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -18,7 +18,17 @@ chapterPreamble() ``` -Building upon the concepts and workflows presented in the previous chapter on mediation analysis, here we demonstrate how to perform **multimodal mediation** **analysis** to identify potential mediators across multiple omics layers. **Multimodal mediation analysis** examines whether and how an exposure (e.g., treatment, diet) affects an outcome (e.g., disease, phenotype) **through** **intermediate variables (mediators) across multiple omics layers** **simultaneously** (e.g., microbiome, metabolomics). It identifies which features across modalities act as mediators, quantifies indirect (mediated) effects while adjusting for covariates, and compares mediation strength across layers, providing mechanistic insights into **how exposures impact outcomes** **via molecular pathways**. +Building upon the concepts and workflows presented in the previous chapter on +mediation analysis, here we demonstrate how to perform **multimodal mediation** +**analysis** to identify potential mediators across multiple omics layers. +**Multimodal mediation analysis** examines whether and how an exposure (e.g., +treatment, diet) affects an outcome (e.g., disease, phenotype) **through** +**intermediate variables (mediators) across multiple omics layers** +**simultaneously** (e.g., microbiome, metabolomics). It identifies which +features across modalities act as mediators, quantifies indirect (mediated) +effects while adjusting for covariates, and compares mediation strength across +layers, providing mechanistic insights into **how exposures impact outcomes** +**via molecular pathways**. ```{r} #| label: fig_multimodal_mediation @@ -54,9 +64,22 @@ digraph multimodal_mediation { ") ``` -In this chapter, we demonstrate multimodal mediation analysis using microbiome data from the iHMP IBDMDB study from the R/Bioconductor package curatedMetagenomicData. [@Lloyd-Price2019]. Unlike conventional mediation analysis, which typically focuses on a single mediator or a small set of mediators, we performed multimodal mediation analysis with R/Bioconductor package \[multimedia\] [@Jiang2025] here, which can handle many potential mediators across multiple data modalities. For example, species-level taxonomic abundances, pathways abundances, and metabolomic profiles can all serve as potential mediators. This method allows us to estimate and rank each mediator’s contribution while accounting for correlations between mediators in a high-dimensional setting. - -In this example, time point (baseline vs. post-treatment) is used as the treatment, and microbiome-based dysbiosis scores are used as the outcome. We will first focus on the relative abundances of gut microbial species as mediators, followed by pathways abundances. +In this chapter, we demonstrate multimodal mediation analysis using microbiome +data from the iHMP IBDMDB study from the R/Bioconductor package +curatedMetagenomicData. [@Lloyd-Price2019]. Unlike conventional mediation +analysis, which typically focuses on a single mediator or a small set of +mediators, we performed multimodal mediation analysis with R/Bioconductor +package \[multimedia\] [@Jiang2025] here, which can handle many potential +mediators across multiple data modalities. For example, species-level taxonomic +abundances, pathways abundances, and metabolomic profiles can all serve as +potential mediators. This method allows us to estimate and rank each mediator’s +contribution while accounting for correlations between mediators in a high- +dimensional setting. + +In this example, time point (baseline vs. post-treatment) is used as the +treatment, and microbiome-based dysbiosis scores are used as the outcome. We +will first focus on the relative abundances of gut microbial species as +mediators, followed by pathways abundances. Generally, we will proceed through the following key steps: @@ -68,15 +91,21 @@ Generally, we will proceed through the following key steps: 4. Fit the multimodal mediation model using the R package multimedia. -5. Interpret both the overall indirect effects and the mediator-specific indirect effects. +5. Interpret both the overall indirect effects and the mediator-specific +indirect effects. -6. Visualize the results using forest plots, histograms, and rankings to prioritize findings. +6. Visualize the results using forest plots, histograms, and rankings to +prioritize findings. 7. Repeat the process for the iHMP microbial pathways. ### Performing multimodal mediation analysis for iHMP species relative abundance {#sec-relative-abundance} -We begin by loading the relative abundance data and calculating dysbiosis scores for each sample based on Bray–Curtis dissimilarity from a healthy reference set. This represents deviation from a healthy microbiome, where higher dysbiosis scores indicate greater microbial imbalance. We then subset to IBD subjects with samples collected at both baseline (visit 1) and post-treatment (visit 21). +We begin by loading the relative abundance data and calculating dysbiosis scores +for each sample based on Bray–Curtis dissimilarity from a healthy reference set. +This represents deviation from a healthy microbiome, where higher dysbiosis +scores indicate greater microbial imbalance. We then subset to IBD subjects with +samples collected at both baseline (visit 1) and post-treatment (visit 21). ```{r} #| label: iHMP_data_preprocessing @@ -85,8 +114,7 @@ We begin by loading the relative abundance data and calculating dysbiosis scores ################## # Load libraries # ################## -remove(list = ls()) -gc() + library(curatedMetagenomicData) library(SummarizedExperiment) library(dplyr) @@ -102,43 +130,31 @@ library(mia) ################## # Load data -tse_relative <- curatedMetagenomicData( +tse <- curatedMetagenomicData( "HMP_2019_ibdmdb.relative_abundance", + rownames = "short", dryrun = FALSE )[[1]] # Assign SampleID for matching -colData(tse_relative)$SampleID <- colnames(tse_relative) +tse[["SampleID"]] <- colnames(tse) # Convert relative_abundance assay to relabundance (which is in [0,1] interval) -tse_relative <- transformAssay(tse_relative, assay.type="relative_abundance", method="relabundance") - -# Optionally, remove the original assay to avoid confusion between the two relative abundance versions -assay(tse_relative, "relative_abundance") <- NULL - -# Remove samples with NA in the relabundance assay -tse_relative <- tse_relative[, colSums(is.na(assay(tse_relative))) == 0] - -# Change the rownames names for variables of interest -rownames(tse_relative) <- sub('.*s__', '', rownames(tse_relative)) -rownames(tse_relative) <- str_remove_all(rownames(tse_relative), "\\[|\\]") -rownames(tse_relative) <- str_replace_all(rownames(tse_relative), "[: \\.,]", "_") -safe_names <- make.names(rownames(tse_relative), unique = TRUE) -rownames(tse_relative) <- safe_names +tse <- transformAssay(tse, assay.type="relative_abundance", method="relabundance") ######################## # Reference nonIBD set # ######################## # Reference set: healthy -tse_relative$disease_binary <- tse_relative$disease == "healthy" +tse$disease_binary <- tse$disease == "healthy" ######################################## # Calculate Bray-Curtis dissimilarity # ######################################## # Bray-Curtis dissimilarity -diss <- as.matrix(getDissimilarity(tse_relative, method = "bray", na.rm=TRUE, assay.type = "relabundance")) +diss <- as.matrix(getDissimilarity(tse, method = "bray", na.rm=TRUE, assay.type = "relabundance")) ################################# # Calculate the dysbiosis score # @@ -146,12 +162,12 @@ diss <- as.matrix(getDissimilarity(tse_relative, method = "bray", na.rm=TRUE, as # Calculate dysbiosis score for each sample # For each sample i, we compute the median distance between sample i and all reference samples in `ref_set` -sample_ids <- colData(tse_relative)$SampleID +sample_ids <- tse[["SampleID"]] -colData(tse_relative)$dysbiosis <- sapply(seq_along(tse_relative$disease_binary), function(i) { +tse[["dysbiosis"]] <- sapply(seq_along(tse$disease_binary), function(i) { # Logical vector indicating all other reference samples (excluding i) - ref_others <- tse_relative$disease_binary & (sample_ids != sample_ids[i]) + ref_others <- tse$disease_binary & (sample_ids != sample_ids[i]) # Compute median distance between sample i and these reference samples median(diss[i, ref_others], na.rm = TRUE) @@ -163,32 +179,36 @@ colData(tse_relative)$dysbiosis <- sapply(seq_along(tse_relative$disease_binary) ############################################################# # Keep IBD only -tse_relative <- tse_relative[, colData(tse_relative)$disease != "healthy"] +tse <- tse[, tse[["disease"]] != "healthy"] # Visit 1 (baseline) or 25 (post) -tse_relative <- tse_relative[, colData(tse_relative)$visit_number %in% c(1, 21)] +tse <- tse[, tse[["visit_number"]] %in% c(1, 21)] # Keep subjects with both visits -keep_subjects <- names(which(table(colData(tse_relative)$subject_id) == 2)) -tse_relative <- tse_relative[, colData(tse_relative)$subject_id %in% keep_subjects] +keep_subjects <- names(which(table(tse[["subject_id"]]) == 2)) +tse <- tse[, tse[["subject_id"]] %in% keep_subjects] # Define time point label -colData(tse_relative)$Time_point <- ifelse( - colData(tse_relative)$visit_number == 1, "T0", "T1" +tse[["Time_point"]] <- ifelse( + tse[["visit_number"]] == 1, "T0", "T1" ) # Define treatment: Pre and Post (Baseline vs. Post-baseline) -colData(tse_relative)$treatment <- factor(colData(tse_relative)$Time_point, levels = c("T0", "T1")) +tse[["treatment"]] <- factor(tse[["Time_point"]], levels = c("T0", "T1")) ``` -Next, we will define the mediation analysis components. By using the processed data, we define: +Next, we will define the mediation analysis components. By using the processed +data, we define: - Treatment: time point (baseline vs. post-baseline) - Outcome: dysbiosis score - Mediators: scaled species abundances -These are bundled into a mediation data object (exper), which is passed into the multimedia function to estimate path models for the treatment-to-mediator (α path), mediator-to-outcome (β path), and the combined indirect effect (α×β path). +These are bundled into a mediation data object (exper), which is passed into the +multimedia function to estimate path models for the treatment-to-mediator (α +path), mediator-to-outcome (β path), and the combined indirect effect (α×β +path). ```{r} #| label: define_mediation_data_frame @@ -198,16 +218,14 @@ These are bundled into a mediation data object (exper), which is passed into the # Define mediation data set # ############################# -# Outcome: Dysbiosis Score -colData(tse_relative)$dysbiosis <- as.numeric(colData(tse_relative)$dysbiosis) - # Mediators (scaled, directly from SE): Species -tse_relative <- transformAssay(tse_relative, assay.type="relabundance", method="standardize", name="scaled") -assays(tse_relative) <- SimpleList(mediators = assay(tse_relative, "scaled")) +tse <- transformAssay(tse, assay.type="relabundance", method="standardize", name="scaled") ``` -We now fit the multimodal mediation model and inspect the overall indirect and direct effects to understand whether microbiome composition mediates the treatment's effect on dysbiosis. +We now fit the multimodal mediation model and inspect the overall indirect and +direct effects to understand whether microbiome composition mediates the +treatment's effect on dysbiosis. ```{r} #| label: overall_mediation_analysis @@ -218,19 +236,15 @@ We now fit the multimodal mediation model and inspect the overall indirect and d ############################# # Convert the TreeSummarizedExperiment to be SummarizedExperiment -se_relative <- tryCatch(as(tse_relative, "SummarizedExperiment"), error = function(e) NULL) -if (is.null(rownames(se_relative)) || any(rownames(se_relative) == "")) { - rn <- rownames(tse_relative) - stopifnot(length(rn) == nrow(se_relative)) - rownames(se_relative) <- rn - rownames(rowData(se_relative)) <- rn # keep rowData aligned too -} +se_relative <- as(tse, "SummarizedExperiment") +rownames(se_relative) <- rownames(tse) +rownames(rowData(se_relative)) <- rownames(tse) # keep rowData aligned too # Use indices for mediators medi_idx <- seq_len(nrow(se_relative)) # Create the Mediation Data object -exper <- multimedia::mediation_data( +exper <- mediation_data( se_relative, outcomes = "dysbiosis", treatments = "treatment", @@ -248,7 +262,6 @@ print(direct_effect(res, exper)) ``` - ```{r} #| label: specific_mediation #| message: false @@ -265,7 +278,7 @@ extract_effects <- function(res, treat_term = "treatmentT1", outcome = "dysbiosi # β: mediator -> outcome (single model with all mediators) beta <- unname(coef(res@outcome@estimates[[outcome]])[names(alpha)]) - tibble::tibble( + tibble( mediator = names(alpha), alpha = as.numeric(alpha), beta = as.numeric(beta), @@ -285,7 +298,9 @@ print(top_mediators) ``` -To quantify uncertainty around the indirect effects, we use non-parametric bootstrap resampling to compute **95% confidence intervals** for both the overall and mediator-specific indirect effects. +To quantify uncertainty around the indirect effects, we use non-parametric +bootstrap resampling to compute **95% confidence intervals** for both the +overall and mediator-specific indirect effects. ```{r} #| label: visualization_for_overall_indirect_effects @@ -312,7 +327,9 @@ ggplot(boot_overall$indirect) + ``` -Finally, we visualize the **mediator-specific indirect effects** using a forest plot, displaying point estimates and 95% CIs to highlight the most influential mediators. +Finally, we visualize the **mediator-specific indirect effects** using a forest +plot, displaying point estimates and 95% CIs to highlight the most influential +mediators. ```{r} #| label: visualization_for_mediation_specific_indirect_effects @@ -378,49 +395,56 @@ ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=signific ``` -Based on the forest plot results, we can see that the Agathobaculum butyriciproducens and Escherichia coli pathways has significant indirect mediation effects. +Based on the forest plot results, we can see that the Agathobaculum +butyriciproducens and Escherichia coli pathways has significant indirect +mediation effects. -Next we repeat the mediation analysis with the iHMP pathway abundances using the same setup as before, using the time point as the treatment, dysbiosis score as the outcome, but pathways abundances as the mediation variables instead of the relative species abundance. +Next we repeat the mediation analysis with the iHMP pathway abundances using the +same setup as before, using the time point as the treatment, dysbiosis score as +the outcome, but pathways abundances as the mediation variables instead of the +relative species abundance. ### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} -Let's extract the pathway abundance data from the iHMP data first, as well as other data pre-processing. +Let's extract the pathway abundance data from the iHMP data first, as well as +other data pre-processing. ```{r} #| label: load_pkg_data #| message: false -tse_pathway <- curatedMetagenomicData( +tse <- curatedMetagenomicData( "HMP_2019_ibdmdb.pathway_abundance", + rownames = "short", dryrun = FALSE )[[1]] -colData(tse_pathway)$SampleID <- colnames(tse_pathway) -tse_pathway <- tse_pathway[, colSums(is.na(assay(tse_pathway))) == 0] +tse[["SampleID"]] <- colnames(tse) +tse <- tse[, colSums(is.na(assay(tse))) == 0] -rows_to_keep <- !grepl("\\|", rownames(tse_pathway)) +rows_to_keep <- !grepl("\\|", rownames(tse)) rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE -tse_pathway <- tse_pathway[rows_to_keep, ] -assay(tse_pathway) <- assay(tse_pathway) / 100 -tse_pathway$disease_binary <- colData(tse_pathway)$disease == "healthy" - -diss <- as.matrix(vegdist(t(assay(tse_pathway)), method = "bray", na.rm = TRUE)) -colData(tse_pathway)$dysbiosis <- sapply(seq_along(tse_pathway$disease_binary), function(i) { - median(diss[i, tse_pathway$disease_binary & - (colData(tse_pathway)$SampleID != colData(tse_pathway)$SampleID[i])], +tse <- tse[rows_to_keep, ] +assay(tse) <- assay(tse) / 100 +tse$disease_binary <- tse[["disease"]] == "healthy" + +diss <- as.matrix(vegdist(t(assay(tse)), method = "bray", na.rm = TRUE)) +tse[["dysbiosis"]] <- sapply(seq_along(tse$disease_binary), function(i) { + median(diss[i, tse$disease_binary & + (tse[["SampleID"]] != tse[["SampleID"]][i])], na.rm = TRUE) }) -tse_pathway <- tse_pathway[, colData(tse_pathway)$disease != "healthy"] -tse_pathway <- tse_pathway[, colData(tse_pathway)$visit_number %in% c(1, 27)] -keep_subjects <- names(which(table(colData(tse_pathway)$subject_id) == 2)) -tse_pathway <- tse_pathway[, colData(tse_pathway)$subject_id %in% keep_subjects] +tse <- tse[, tse[["disease"]] != "healthy"] +tse <- tse[, tse[["visit_number"]] %in% c(1, 27)] +keep_subjects <- names(which(table(tse[["subject_id"]]) == 2)) +tse <- tse[, tse[["subject_id"]] %in% keep_subjects] -colData(tse_pathway)$Time_point <- ifelse( - colData(tse_pathway)$visit_number == 1, "T0", "T1" +tse[["Time_point"]] <- ifelse( + tse[["visit_number"]] == 1, "T0", "T1" ) -colData(tse_pathway)$treatment <- factor(colData(tse_pathway)$Time_point, levels = c("T0", "T1")) +tse[["treatment"]] <- factor(tse[["Time_point"]], levels = c("T0", "T1")) ``` @@ -430,25 +454,21 @@ We will define the mediation analysis data set then. #| label: define_mediation_data_frame #| message: false -colData(tse_pathway)$dysbiosis <- as.numeric(colData(tse_pathway)$dysbiosis) -tse_pathway <- transformAssay(tse_pathway, assay.type="pathway_abundance", method="standardize", name="scaled") -assays(tse_pathway) <- SimpleList(mediators = assay(tse_pathway, "scaled")) +tse <- transformAssay(tse, assay.type="pathway_abundance", method="standardize", name="scaled") ``` -Next, we continue to fit the multimodal mediation model, and extract the overall and mediation-specific indirect effects. +Next, we continue to fit the multimodal mediation model, and extract the overall +and mediation-specific indirect effects. ```{r} #| label: mediation_analysis #| message: false -se_pathway <- tryCatch(as(tse_pathway, "SummarizedExperiment"), error = function(e) NULL) -if (is.null(rownames(se_pathway)) || any(rownames(se_pathway) == "")) { - rn <- rownames(tse_pathway) - stopifnot(length(rn) == nrow(se_pathway)) - rownames(se_pathway) <- rn - rownames(rowData(se_pathway)) <- rn # keep rowData aligned too -} +se_pathway <- as(tse, "SummarizedExperiment") +rownames(se_pathway) <- rownames(tse) +rownames(rowData(se_pathway)) <- rownames(tse) # keep rowData aligned too + raw <- rownames(se_pathway) clean1 <- str_remove_all(raw, "\\[|\\]") @@ -458,7 +478,7 @@ rownames(se_pathway) <- safe_names medi_idx <- seq_len(nrow(se_pathway)) -exper <- multimedia::mediation_data( +exper <- mediation_data( se_pathway, outcomes = "dysbiosis", treatments = "treatment", @@ -480,7 +500,8 @@ print(top_mediators) ``` -The non-parametric bootstrap resampling is also used to compute **95% CIs** for both the overall and mediator-specific pathway indirect effects. +The non-parametric bootstrap resampling is also used to compute **95% CIs** for +both the overall and mediator-specific pathway indirect effects. ```{r} #| label: visualization_for_overall_indirect_effects @@ -502,7 +523,8 @@ ggplot(boot_overall$indirect) + ``` -We visualize the **mediator-specific indirect effects** using a forest plot for pathway aubundance as well. with the point estimates and 95% CIs. +We visualize the **mediator-specific indirect effects** using a forest plot for +pathway aubundance as well. with the point estimates and 95% CIs. ```{r} #| label: visualization_for_mediation_specific_indirect_effects @@ -563,4 +585,5 @@ ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=signific ``` -Based on the forest plot results, we can see that the tRNA charging pathway has significant indirect mediation effects on the dysbiosis scores. +Based on the forest plot results, we can see that the tRNA charging pathway has +significant indirect mediation effects on the dysbiosis scores. From 00fe5faf067f5bbca1ec3136a00e677190ee504e Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Mon, 29 Dec 2025 20:59:30 +0800 Subject: [PATCH 06/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 122 ++++++++++++++++++-------------------- 1 file changed, 59 insertions(+), 63 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 18ef63f1..f5b9bd8f 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -137,7 +137,7 @@ tse <- curatedMetagenomicData( )[[1]] # Assign SampleID for matching -tse[["SampleID"]] <- colnames(tse) +tse$SampleID <- colnames(tse) # Convert relative_abundance assay to relabundance (which is in [0,1] interval) tse <- transformAssay(tse, assay.type="relative_abundance", method="relabundance") @@ -161,17 +161,15 @@ diss <- as.matrix(getDissimilarity(tse, method = "bray", na.rm=TRUE, assay.type ################################# # Calculate dysbiosis score for each sample -# For each sample i, we compute the median distance between sample i and all reference samples in `ref_set` -sample_ids <- tse[["SampleID"]] - -tse[["dysbiosis"]] <- sapply(seq_along(tse$disease_binary), function(i) { +# For each sample i, we compute the median distance between sample i and all reference samples +sample_ids <- tse$SampleID +tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { # Logical vector indicating all other reference samples (excluding i) ref_others <- tse$disease_binary & (sample_ids != sample_ids[i]) # Compute median distance between sample i and these reference samples median(diss[i, ref_others], na.rm = TRUE) - }) ############################################################# @@ -179,22 +177,20 @@ tse[["dysbiosis"]] <- sapply(seq_along(tse$disease_binary), function(i) { ############################################################# # Keep IBD only -tse <- tse[, tse[["disease"]] != "healthy"] +tse <- tse[, tse$disease != "healthy"] -# Visit 1 (baseline) or 25 (post) -tse <- tse[, tse[["visit_number"]] %in% c(1, 21)] +# Visit 1 (baseline) or 21 (post) +tse <- tse[, tse$visit_number %in% c(1, 21)] # Keep subjects with both visits -keep_subjects <- names(which(table(tse[["subject_id"]]) == 2)) -tse <- tse[, tse[["subject_id"]] %in% keep_subjects] +keep_subjects <- names(which(table(tse$subject_id) == 2)) +tse <- tse[, tse$subject_id %in% keep_subjects] # Define time point label -tse[["Time_point"]] <- ifelse( - tse[["visit_number"]] == 1, "T0", "T1" -) +tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") # Define treatment: Pre and Post (Baseline vs. Post-baseline) -tse[["treatment"]] <- factor(tse[["Time_point"]], levels = c("T0", "T1")) +tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) ``` @@ -235,17 +231,12 @@ treatment's effect on dysbiosis. # Run mediation (non-delta) # ############################# -# Convert the TreeSummarizedExperiment to be SummarizedExperiment -se_relative <- as(tse, "SummarizedExperiment") -rownames(se_relative) <- rownames(tse) -rownames(rowData(se_relative)) <- rownames(tse) # keep rowData aligned too - # Use indices for mediators -medi_idx <- seq_len(nrow(se_relative)) +medi_idx <- seq_len(nrow(tse)) # Create the Mediation Data object exper <- mediation_data( - se_relative, + tse, outcomes = "dysbiosis", treatments = "treatment", mediators = medi_idx @@ -270,29 +261,13 @@ print(direct_effect(res, exper)) # Mediator-specific effects ############################# -# Helper function to extract the effects -extract_effects <- function(res, treat_term = "treatmentT1", outcome = "dysbiosis") { - # α: treatment -> each mediator (one model per mediator) - alpha <- sapply(res@mediation@estimates, function(m) unname(coef(m)[treat_term])) - - # β: mediator -> outcome (single model with all mediators) - beta <- unname(coef(res@outcome@estimates[[outcome]])[names(alpha)]) - - tibble( - mediator = names(alpha), - alpha = as.numeric(alpha), - beta = as.numeric(beta), - indirect_effect = alpha * beta - ) -} - # Extract the treatment-to-mediator path coefficients -effects_by_mediator <- extract_effects(res, treat_term = "treatmentT1", outcome = "dysbiosis") +effects_by_mediator <- indirect_pathwise(res, exper) # Rank top_mediators <- effects_by_mediator %>% arrange(desc(abs(indirect_effect))) %>% - head(10) + head(20) print(top_mediators) @@ -404,6 +379,33 @@ same setup as before, using the time point as the treatment, dysbiosis score as the outcome, but pathways abundances as the mediation variables instead of the relative species abundance. +```{r} +#| label: default_visualization_plot_mediators_pathway +#| message: false + +# pick top mediators from the already-computed bootstrap summary (fast) +top_meds <- summary_df %>% + dplyr::arrange(dplyr::desc(abs(estimate))) %>% + dplyr::slice_head(n = 12) %>% + dplyr::pull(mediator) + +# construct a minimal effect table that plot_mediators() can use +ie_pw_fast <- tibble::tibble( + outcome = "dysbiosis", + mediator = top_meds, + direct_setting = levels(tse$treatment)[1], + contrast = paste(levels(tse$treatment)[1], "-", levels(tse$treatment)[2]), + indirect_effect = summary_df$estimate[match(top_meds, summary_df$mediator)] +) + +plot_mediators(ie_pw_fast, exper, n_panels = 12) + + +``` +Mediators are standardized (z-scored) for comparability across features; hence values may be negative. Besides, many microbial features are sparse; several mediators exhibit near-zero values for most samples, producing vertical bands. + + + ### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} Let's extract the pathway abundance data from the iHMP data first, as well as @@ -419,32 +421,30 @@ tse <- curatedMetagenomicData( dryrun = FALSE )[[1]] -tse[["SampleID"]] <- colnames(tse) +tse$SampleID <- colnames(tse) tse <- tse[, colSums(is.na(assay(tse))) == 0] rows_to_keep <- !grepl("\\|", rownames(tse)) rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE tse <- tse[rows_to_keep, ] assay(tse) <- assay(tse) / 100 -tse$disease_binary <- tse[["disease"]] == "healthy" +tse$disease_binary <- tse$disease == "healthy" diss <- as.matrix(vegdist(t(assay(tse)), method = "bray", na.rm = TRUE)) -tse[["dysbiosis"]] <- sapply(seq_along(tse$disease_binary), function(i) { +tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { median(diss[i, tse$disease_binary & - (tse[["SampleID"]] != tse[["SampleID"]][i])], + (tse$SampleID != tse$SampleID[i])], na.rm = TRUE) }) -tse <- tse[, tse[["disease"]] != "healthy"] -tse <- tse[, tse[["visit_number"]] %in% c(1, 27)] -keep_subjects <- names(which(table(tse[["subject_id"]]) == 2)) -tse <- tse[, tse[["subject_id"]] %in% keep_subjects] +tse <- tse[, tse$disease != "healthy"] +tse <- tse[, tse$visit_number %in% c(1, 27)] +keep_subjects <- names(which(table(tse$subject_id) == 2)) +tse <- tse[, tse$subject_id %in% keep_subjects] -tse[["Time_point"]] <- ifelse( - tse[["visit_number"]] == 1, "T0", "T1" -) +tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") -tse[["treatment"]] <- factor(tse[["Time_point"]], levels = c("T0", "T1")) +tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) ``` @@ -465,21 +465,16 @@ and mediation-specific indirect effects. #| label: mediation_analysis #| message: false -se_pathway <- as(tse, "SummarizedExperiment") -rownames(se_pathway) <- rownames(tse) -rownames(rowData(se_pathway)) <- rownames(tse) # keep rowData aligned too - - -raw <- rownames(se_pathway) +raw <- rownames(tse) clean1 <- str_remove_all(raw, "\\[|\\]") clean2 <- str_replace_all(clean1, "[: \\.,]", "_") safe_names <- make.names(clean2, unique = TRUE) -rownames(se_pathway) <- safe_names +rownames(tse) <- safe_names -medi_idx <- seq_len(nrow(se_pathway)) +medi_idx <- seq_len(nrow(tse)) exper <- mediation_data( - se_pathway, + tse, outcomes = "dysbiosis", treatments = "treatment", mediators = medi_idx @@ -491,11 +486,12 @@ summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) -effects_by_mediator <- extract_effects(res, treat_term = "treatmentT1", outcome = "dysbiosis") +effects_by_mediator <- indirect_pathwise(res, exper) top_mediators <- effects_by_mediator %>% arrange(desc(abs(indirect_effect))) %>% - head(10) + head(20) + print(top_mediators) ``` From 2c1b06cb633d33988700a6539686eef01517532c Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Mon, 26 Jan 2026 07:22:25 +0800 Subject: [PATCH 07/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index f5b9bd8f..a16d5118 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -165,11 +165,13 @@ diss <- as.matrix(getDissimilarity(tse, method = "bray", na.rm=TRUE, assay.type sample_ids <- tse$SampleID tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { + # Logical vector indicating all other reference samples (excluding i) ref_others <- tse$disease_binary & (sample_ids != sample_ids[i]) # Compute median distance between sample i and these reference samples median(diss[i, ref_others], na.rm = TRUE) + }) ############################################################# @@ -231,12 +233,16 @@ treatment's effect on dysbiosis. # Run mediation (non-delta) # ############################# +# Convert the TreeSummarizedExperiment to be SummarizedExperiment +se_relative <- as(tse, "SummarizedExperiment") +if (is.null(rownames(se_relative))) rownames(se_relative) <- rownames(rowData(se_relative)) <- make.names(rownames(tse), unique = TRUE) + # Use indices for mediators -medi_idx <- seq_len(nrow(tse)) +medi_idx <- seq_len(nrow(se_relative)) # Create the Mediation Data object exper <- mediation_data( - tse, + se_relative, outcomes = "dysbiosis", treatments = "treatment", mediators = medi_idx @@ -286,7 +292,7 @@ overall and mediator-specific indirect effects. ##################################### set.seed(12345) -boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) +boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 10) summary(boot_overall$indirect) quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) From eea2ec11a0146b998e0da28cf937ad135d374c2c Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Sun, 19 Jul 2026 01:12:37 +0800 Subject: [PATCH 08/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 443 ++++++++++++++------------------------ 1 file changed, 159 insertions(+), 284 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index a16d5118..711f57a2 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -40,9 +40,9 @@ library(DiagrammeR) grViz(" digraph multimodal_mediation { graph [layout = dot, rankdir = LR] - + node [fontname = Helvetica, fontsize = 12] - + Exposure [shape=box, style=filled, fillcolor=lightblue, color=black] Outcome [shape=box, style=filled, fillcolor=palegreen, color=black] @@ -107,14 +107,12 @@ This represents deviation from a healthy microbiome, where higher dysbiosis scores indicate greater microbial imbalance. We then subset to IBD subjects with samples collected at both baseline (visit 1) and post-treatment (visit 21). +We first load the packages and define small plotting helpers used below. + ```{r} -#| label: iHMP_data_preprocessing +#| label: load_libraries_and_helpers #| message: false -################## -# Load libraries # -################## - library(curatedMetagenomicData) library(SummarizedExperiment) library(dplyr) @@ -124,80 +122,74 @@ library(multimedia) library(stringr) library(ggplot2) library(mia) - -################## -# Load iHMP data # -################## - -# Load data -tse <- curatedMetagenomicData( - "HMP_2019_ibdmdb.relative_abundance", - rownames = "short", - dryrun = FALSE - )[[1]] - -# Assign SampleID for matching -tse$SampleID <- colnames(tse) - -# Convert relative_abundance assay to relabundance (which is in [0,1] interval) -tse <- transformAssay(tse, assay.type="relative_abundance", method="relabundance") - -######################## -# Reference nonIBD set # -######################## - -# Reference set: healthy -tse$disease_binary <- tse$disease == "healthy" - -######################################## -# Calculate Bray-Curtis dissimilarity # -######################################## - -# Bray-Curtis dissimilarity -diss <- as.matrix(getDissimilarity(tse, method = "bray", na.rm=TRUE, assay.type = "relabundance")) - -################################# -# Calculate the dysbiosis score # -################################# - -# Calculate dysbiosis score for each sample -# For each sample i, we compute the median distance between sample i and all reference samples -sample_ids <- tse$SampleID - -tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { - - # Logical vector indicating all other reference samples (excluding i) - ref_others <- tse$disease_binary & (sample_ids != sample_ids[i]) +library(miaViz) + +plot_indirect_histogram <- function(boot_indirect, title) { + histogram_data <- SummarizedExperiment( + assays = list( + indirect_effect = matrix( + boot_indirect$indirect_effect, + nrow = 1, + dimnames = list("indirect_effect", seq_len(nrow(boot_indirect))) + ) + ) + ) - # Compute median distance between sample i and these reference samples - median(diss[i, ref_others], na.rm = TRUE) - -}) + miaViz::plotHistogram(histogram_data, assay.type = "indirect_effect") + + labs( + x = "Overall Indirect Effect", + y = "Frequency", + title = title + ) +} -############################################################# -# Subset to IBD only and only one post-baseline time point # -############################################################# +plot_indirect_forest <- function(summary_df, title) { + forest_data <- summary_df %>% + mutate(abs_estimate = abs(estimate)) + + miaViz::plotForest( + forest_data, + effect.var = "estimate", + ci.lower.var = "lower", + ci.upper.var = "upper", + pval.var = NULL, + id.var = "mediator", + order.by = "abs_estimate", + colour.by = "significant" + ) + + labs( + x = "Observed Indirect Effect with Bootstrap CI", + y = "Mediator", + title = title + ) +} -# Keep IBD only -tse <- tse[, tse$disease != "healthy"] +clean_pathway_names <- function(pathways) { + pathways <- gsub("PWY0\\.", "PWY0-", pathways) + pathways <- gsub("PWY\\.", "PWY-", pathways) + pathways <- gsub("__", ": ", pathways) + pathways <- gsub("_\\.", " ", pathways) + pathways <- gsub("\\.", " ", pathways) + pathways <- gsub("_", " ", pathways) + pathways <- gsub("\\s+", " ", pathways) + trimws(pathways) +} -# Visit 1 (baseline) or 21 (post) -tse <- tse[, tse$visit_number %in% c(1, 21)] +``` -# Keep subjects with both visits -keep_subjects <- names(which(table(tse$subject_id) == 2)) -tse <- tse[, tse$subject_id %in% keep_subjects] +The processed demo data are loaded directly from the local data file prepared +for this chapter. -# Define time point label -tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") +```{r} +#| label: load_iHMP_relative_demo_data +#| message: false -# Define treatment: Pre and Post (Baseline vs. Post-baseline) -tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) +load("/Users/lyh_ciel/Library/CloudStorage/OneDrive-YaleUniversity/WCM/OMA/multimedia_data_demo.rda") +tse <- tse_relative ``` -Next, we will define the mediation analysis components. By using the processed -data, we define: +The mediation analysis components are defined as following: - Treatment: time point (baseline vs. post-baseline) - Outcome: dysbiosis score @@ -208,19 +200,6 @@ multimedia function to estimate path models for the treatment-to-mediator (α path), mediator-to-outcome (β path), and the combined indirect effect (α×β path). -```{r} -#| label: define_mediation_data_frame -#| message: false - -############################# -# Define mediation data set # -############################# - -# Mediators (scaled, directly from SE): Species -tse <- transformAssay(tse, assay.type="relabundance", method="standardize", name="scaled") - -``` - We now fit the multimodal mediation model and inspect the overall indirect and direct effects to understand whether microbiome composition mediates the treatment's effect on dysbiosis. @@ -229,51 +208,36 @@ treatment's effect on dysbiosis. #| label: overall_mediation_analysis #| message: false -############################# -# Run mediation (non-delta) # -############################# - -# Convert the TreeSummarizedExperiment to be SummarizedExperiment -se_relative <- as(tse, "SummarizedExperiment") -if (is.null(rownames(se_relative))) rownames(se_relative) <- rownames(rowData(se_relative)) <- make.names(rownames(tse), unique = TRUE) - -# Use indices for mediators medi_idx <- seq_len(nrow(se_relative)) -# Create the Mediation Data object exper <- mediation_data( - se_relative, - outcomes = "dysbiosis", - treatments = "treatment", - mediators = medi_idx + se_relative, + outcomes = "dysbiosis", + treatments = "treatment", + mediators = medi_idx ) -# Fit and summarize mdl <- multimedia(exper) res <- estimate(mdl, exper) -# Summarize overall indirect/direct summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) ``` +Mediator-specific indirect effects are then extracted and ranked by absolute +effect size. + ```{r} #| label: specific_mediation #| message: false -############################# -# Mediator-specific effects -############################# - -# Extract the treatment-to-mediator path coefficients effects_by_mediator <- indirect_pathwise(res, exper) -# Rank top_mediators <- effects_by_mediator %>% - arrange(desc(abs(indirect_effect))) %>% - head(20) + arrange(desc(abs(indirect_effect))) %>% + head(20) print(top_mediators) @@ -284,27 +248,18 @@ bootstrap resampling to compute **95% confidence intervals** for both the overall and mediator-specific indirect effects. ```{r} -#| label: visualization_for_overall_indirect_effects +#| label: visualization_for_species_overall_indirect_effects #| message: false -##################################### -# Bootstrap overall indirect effect # -##################################### - set.seed(12345) boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 10) summary(boot_overall$indirect) quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) -# Histogram -ggplot(boot_overall$indirect) + - geom_histogram(aes(indirect_effect), bins = 20, fill = "#69b3a2", color = "black") + - theme_classic() + - labs( - x = "Overall Indirect Effect", - y = "Frequency", - title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" - ) +plot_indirect_histogram( + boot_overall$indirect, + "Bootstrap Distribution of Overall Indirect Effect (B = 100)" +) ``` @@ -313,66 +268,49 @@ plot, displaying point estimates and 95% CIs to highlight the most influential mediators. ```{r} -#| label: visualization_for_mediation_specific_indirect_effects +#| label: visualization_for_species_specific_indirect_effects #| message: false -############################################### -# Bootstrap mediator-specific indirect effect # -############################################### - -# Define an mediation-specific indirect effect function manually indirect_each <- function(mdl, exper) { - res <- estimate(mdl, exper) - alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) - names(alpha) <- sub("\\.treatmentT1$", "", names(alpha)) - beta <- coef(res@outcome@estimates$dysbiosis)[names(alpha)] - indirect <- alpha * beta - names(indirect) <- names(alpha) - indirect <- indirect[!is.na(indirect)] - return(indirect) + res <- estimate(mdl, exper) + alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) + names(alpha) <- sub("\\.treatmentT1$", "", names(alpha)) + beta <- coef(res@outcome@estimates$dysbiosis)[names(alpha)] + indirect <- alpha * beta + names(indirect) <- names(alpha) + indirect <- indirect[!is.na(indirect)] + return(indirect) } set.seed(12345) boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) summary(boot_each$indirect) -# Summarize the bootstrap results -lower_upper <- apply(boot_each$indirect, 2, function(x) quantile(x, c(0.025, 0.975), na.rm=TRUE)) -means <- apply(boot_each$indirect, 2, mean, na.rm=TRUE) +lower_upper <- apply( + boot_each$indirect, + 2, + function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE) +) +means <- apply(boot_each$indirect, 2, mean, na.rm = TRUE) -# Convert to data frame summary_df <- data.frame( - mediator = colnames(boot_each$indirect), - estimate = means, - lower = lower_upper[1,], - upper = lower_upper[2,] + mediator = colnames(boot_each$indirect), + estimate = means, + lower = lower_upper[1, ], + upper = lower_upper[2, ] ) summary_df <- summary_df[-1, ] -# Mark significance if CI does not cross zero summary_df$significant <- with(summary_df, lower > 0 | upper < 0) - -# Remove the species with 0 confidence interval length -summary_df <- summary_df %>% - filter((upper - lower) != 0) - -# Add rankings -summary_df$rank <- ifelse(summary_df$significant, 1, 2) # 1 for significant, 2 for not +summary_df <- summary_df %>% + filter((upper - lower) != 0) +summary_df$rank <- ifelse(summary_df$significant, 1, 2) summary_df <- summary_df[order(summary_df$rank, -abs(summary_df$estimate)), ] -# Forest plot -ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=significant)) + - geom_point() + - geom_errorbarh(aes(xmin=lower, xmax=upper), height=0.2) + - geom_vline(xintercept=0, linetype="dashed", color="grey") + - theme_classic() + - labs( - x = "Observed Indirect Effect with Bootstrap CI", - y = "Mediator", - title = "Forest Plot of Mediation-specific Indirect Effects for Species" - ) + - scale_color_manual(values=c("black","red")) + - theme(legend.position="bottom") +plot_indirect_forest( + summary_df, + "Forest Plot of Mediation-specific Indirect Effects for Species" +) ``` @@ -380,28 +318,24 @@ Based on the forest plot results, we can see that the Agathobaculum butyriciproducens and Escherichia coli pathways has significant indirect mediation effects. -Next we repeat the mediation analysis with the iHMP pathway abundances using the -same setup as before, using the time point as the treatment, dysbiosis score as -the outcome, but pathways abundances as the mediation variables instead of the -relative species abundance. +We also inspect the top species mediators with the default mediator-level +visualization. ```{r} -#| label: default_visualization_plot_mediators_pathway +#| label: default_visualization_plot_mediators_species #| message: false -# pick top mediators from the already-computed bootstrap summary (fast) top_meds <- summary_df %>% - dplyr::arrange(dplyr::desc(abs(estimate))) %>% - dplyr::slice_head(n = 12) %>% - dplyr::pull(mediator) + dplyr::arrange(dplyr::desc(abs(estimate))) %>% + dplyr::slice_head(n = 12) %>% + dplyr::pull(mediator) -# construct a minimal effect table that plot_mediators() can use ie_pw_fast <- tibble::tibble( - outcome = "dysbiosis", - mediator = top_meds, - direct_setting = levels(tse$treatment)[1], - contrast = paste(levels(tse$treatment)[1], "-", levels(tse$treatment)[2]), - indirect_effect = summary_df$estimate[match(top_meds, summary_df$mediator)] + outcome = "dysbiosis", + mediator = top_meds, + direct_setting = levels(tse$treatment)[1], + contrast = paste(levels(tse$treatment)[1], "-", levels(tse$treatment)[2]), + indirect_effect = summary_df$estimate[match(top_meds, summary_df$mediator)] ) plot_mediators(ie_pw_fast, exper, n_panels = 12) @@ -410,57 +344,21 @@ plot_mediators(ie_pw_fast, exper, n_panels = 12) ``` Mediators are standardized (z-scored) for comparability across features; hence values may be negative. Besides, many microbial features are sparse; several mediators exhibit near-zero values for most samples, producing vertical bands. - +Next we repeat the mediation analysis with the iHMP pathway abundances using the +same setup as before, using the time point as the treatment, dysbiosis score as +the outcome, but pathways abundances as the mediation variables instead of the +relative species abundance. ### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} -Let's extract the pathway abundance data from the iHMP data first, as well as -other data pre-processing. +Let's try to use the pathway abundance data from the iHMP data to do the +analysis as well. ```{r} -#| label: load_pkg_data +#| label: load_iHMP_pathway_demo_data #| message: false -tse <- curatedMetagenomicData( - "HMP_2019_ibdmdb.pathway_abundance", - rownames = "short", - dryrun = FALSE - )[[1]] - -tse$SampleID <- colnames(tse) -tse <- tse[, colSums(is.na(assay(tse))) == 0] - -rows_to_keep <- !grepl("\\|", rownames(tse)) -rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE -tse <- tse[rows_to_keep, ] -assay(tse) <- assay(tse) / 100 -tse$disease_binary <- tse$disease == "healthy" - -diss <- as.matrix(vegdist(t(assay(tse)), method = "bray", na.rm = TRUE)) -tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { - median(diss[i, tse$disease_binary & - (tse$SampleID != tse$SampleID[i])], - na.rm = TRUE) -}) - -tse <- tse[, tse$disease != "healthy"] -tse <- tse[, tse$visit_number %in% c(1, 27)] -keep_subjects <- names(which(table(tse$subject_id) == 2)) -tse <- tse[, tse$subject_id %in% keep_subjects] - -tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") - -tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) - -``` - -We will define the mediation analysis data set then. - -```{r} -#| label: define_mediation_data_frame -#| message: false - -tse <- transformAssay(tse, assay.type="pathway_abundance", method="standardize", name="scaled") +tse <- tse_pathway ``` @@ -471,19 +369,13 @@ and mediation-specific indirect effects. #| label: mediation_analysis #| message: false -raw <- rownames(tse) -clean1 <- str_remove_all(raw, "\\[|\\]") -clean2 <- str_replace_all(clean1, "[: \\.,]", "_") -safe_names <- make.names(clean2, unique = TRUE) -rownames(tse) <- safe_names - medi_idx <- seq_len(nrow(tse)) exper <- mediation_data( - tse, - outcomes = "dysbiosis", - treatments = "treatment", - mediators = medi_idx + tse, + outcomes = "dysbiosis", + treatments = "treatment", + mediators = medi_idx ) mdl <- multimedia(exper) @@ -495,8 +387,8 @@ print(direct_effect(res, exper)) effects_by_mediator <- indirect_pathwise(res, exper) top_mediators <- effects_by_mediator %>% - arrange(desc(abs(indirect_effect))) %>% - head(20) + arrange(desc(abs(indirect_effect))) %>% + head(20) print(top_mediators) @@ -506,7 +398,7 @@ The non-parametric bootstrap resampling is also used to compute **95% CIs** for both the overall and mediator-specific pathway indirect effects. ```{r} -#| label: visualization_for_overall_indirect_effects +#| label: visualization_for_pathway_overall_indirect_effects #| message: false set.seed(1234) @@ -514,76 +406,59 @@ boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) summary(boot_overall$indirect) quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) -ggplot(boot_overall$indirect) + - geom_histogram(aes(indirect_effect), bins = 20, fill = "#69b3a2", color = "black") + - theme_classic() + - labs( - x = "Overall Indirect Effect", - y = "Frequency", - title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" - ) +plot_indirect_histogram( + boot_overall$indirect, + "Bootstrap Distribution of Overall Indirect Effect (B = 100)" +) ``` We visualize the **mediator-specific indirect effects** using a forest plot for -pathway aubundance as well. with the point estimates and 95% CIs. +pathway abundance as well, with the point estimates and 95% CIs. ```{r} -#| label: visualization_for_mediation_specific_indirect_effects +#| label: visualization_for_pathway_specific_indirect_effects #| message: false indirect_each <- function(mdl, exper) { - res <- estimate(mdl, exper) - alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) - names(alpha) <- sub("\\.treatmentT1$", "", names(alpha)) - beta <- coef(res@outcome@estimates$dysbiosis)[names(alpha)] - indirect <- alpha * beta - names(indirect) <- names(alpha) - indirect <- indirect[!is.na(indirect)] - return(indirect) + res <- estimate(mdl, exper) + alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) + names(alpha) <- sub("\\.treatmentT1$", "", names(alpha)) + beta <- coef(res@outcome@estimates$dysbiosis)[names(alpha)] + indirect <- alpha * beta + names(indirect) <- names(alpha) + indirect <- indirect[!is.na(indirect)] + return(indirect) } set.seed(1234) boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) -lower_upper <- apply(boot_each$indirect, 2, function(x) quantile(x, c(0.025, 0.975), na.rm=TRUE)) -means <- apply(boot_each$indirect, 2, mean, na.rm=TRUE) +lower_upper <- apply( + boot_each$indirect, + 2, + function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE) +) +means <- apply(boot_each$indirect, 2, mean, na.rm = TRUE) summary_df <- data.frame( - mediator = colnames(boot_each$indirect), - estimate = means, - lower = lower_upper[1,], - upper = lower_upper[2,] + mediator = colnames(boot_each$indirect), + estimate = means, + lower = lower_upper[1, ], + upper = lower_upper[2, ] ) summary_df <- summary_df[-1, ] summary_df$significant <- with(summary_df, lower > 0 | upper < 0) summary_df <- summary_df[summary_df$upper != summary_df$lower, ] -# Clean the pathway names -pwy <- rownames(summary_df) -pwy <- gsub("PWY0\\.", "PWY0-", pwy) # PWY0.xxx → PWY0-xxx -pwy <- gsub("PWY\\.", "PWY-", pwy) # PWY.xxx → PWY-xxx -pwy <- gsub("\\.", " ", pwy) # leftover dots → spaces -pwy <- gsub("__", ": ", pwy) # double underscores → colon -pwy <- gsub("_\\.", " ", pwy) # underscore then dot → space -pwy <- gsub("_", " ", pwy) # remaining underscores → space -pwy <- gsub("\\.\\.", " ", pwy) # double dots → space -pwy <- gsub("\\.$", "", pwy) # remove ending period -# pwy <- trimws(pwy) # final cleanup +pwy <- clean_pathway_names(rownames(summary_df)) rownames(summary_df) <- pwy +summary_df$mediator <- pwy -ggplot(summary_df, aes(x=estimate, y=reorder(mediator, estimate), color=significant)) + - geom_point() + - geom_errorbarh(aes(xmin=lower, xmax=upper), height=0.2) + - geom_vline(xintercept=0, linetype="dashed", color="grey") + - theme_classic() + - labs( - x = "Observed Indirect Effect with Bootstrap CI", - y = "Mediator", - title = "Forest Plot of Mediation-specific Indirect Effects for Pathway" - ) + - scale_color_manual(values=c("black","red")) + - theme(legend.position="bottom") +plot_indirect_forest( + summary_df, + "Forest Plot of Mediation-specific Indirect Effects for Pathway" +) ``` From 410f5a22b0c993c254425c6266190988780d5cc5 Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Tue, 21 Jul 2026 02:03:53 +0800 Subject: [PATCH 09/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 290 ++++++++++++++++++++++++-------------- 1 file changed, 188 insertions(+), 102 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 711f57a2..f6bad9ab 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -1,10 +1,3 @@ ---- -title: "Multimodal Mediation Analysis" -format: html -editor: visual -author: Yihan Liu, Himel Mallick ---- - # **Multimodal Mediation Analysis** {#sec-MSEA} ```{r} @@ -15,7 +8,6 @@ remove(list = ls()) invisible(gc()) library(rebook) chapterPreamble() - ``` Building upon the concepts and workflows presented in the previous chapter on @@ -39,27 +31,27 @@ library(DiagrammeR) grViz(" digraph multimodal_mediation { - graph [layout = dot, rankdir = LR] - - node [fontname = Helvetica, fontsize = 12] - - Exposure [shape=box, style=filled, fillcolor=lightblue, color=black] - Outcome [shape=box, style=filled, fillcolor=palegreen, color=black] - - M1 [label='Species 1', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] - M2 [label='Species 2', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] - M3 [label='Metabolite 1', shape=ellipse, style=filled, fillcolor=lightpink, color=black] - M4 [label='Metabolite 2', shape=ellipse, style=filled, fillcolor=lightpink, color=black] - - Exposure -> M1 - Exposure -> M2 - Exposure -> M3 - Exposure -> M4 - M1 -> Outcome - M2 -> Outcome - M3 -> Outcome - M4 -> Outcome - Exposure -> Outcome [style=dashed] + graph [layout = dot, rankdir = LR] + + node [fontname = Helvetica, fontsize = 12] + + Exposure [shape=box, style=filled, fillcolor=lightblue, color=black] + Outcome [shape=box, style=filled, fillcolor=palegreen, color=black] + + M1 [label='Species 1', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] + M2 [label='Species 2', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] + M3 [label='Metabolite 1', shape=ellipse, style=filled, fillcolor=lightpink, color=black] + M4 [label='Metabolite 2', shape=ellipse, style=filled, fillcolor=lightpink, color=black] + + Exposure -> M1 + Exposure -> M2 + Exposure -> M3 + Exposure -> M4 + M1 -> Outcome + M2 -> Outcome + M3 -> Outcome + M4 -> Outcome + Exposure -> Outcome [style=dashed] } ") ``` @@ -124,46 +116,6 @@ library(ggplot2) library(mia) library(miaViz) -plot_indirect_histogram <- function(boot_indirect, title) { - histogram_data <- SummarizedExperiment( - assays = list( - indirect_effect = matrix( - boot_indirect$indirect_effect, - nrow = 1, - dimnames = list("indirect_effect", seq_len(nrow(boot_indirect))) - ) - ) - ) - - miaViz::plotHistogram(histogram_data, assay.type = "indirect_effect") + - labs( - x = "Overall Indirect Effect", - y = "Frequency", - title = title - ) -} - -plot_indirect_forest <- function(summary_df, title) { - forest_data <- summary_df %>% - mutate(abs_estimate = abs(estimate)) - - miaViz::plotForest( - forest_data, - effect.var = "estimate", - ci.lower.var = "lower", - ci.upper.var = "upper", - pval.var = NULL, - id.var = "mediator", - order.by = "abs_estimate", - colour.by = "significant" - ) + - labs( - x = "Observed Indirect Effect with Bootstrap CI", - y = "Mediator", - title = title - ) -} - clean_pathway_names <- function(pathways) { pathways <- gsub("PWY0\\.", "PWY0-", pathways) pathways <- gsub("PWY\\.", "PWY-", pathways) @@ -174,19 +126,70 @@ clean_pathway_names <- function(pathways) { pathways <- gsub("\\s+", " ", pathways) trimws(pathways) } - ``` -The processed demo data are loaded directly from the local data file prepared -for this chapter. +We first load the relative abundance data, calculate dysbiosis scores from the +healthy reference samples, and retain IBD subjects with both baseline and +post-treatment samples. ```{r} #| label: load_iHMP_relative_demo_data #| message: false -load("/Users/lyh_ciel/Library/CloudStorage/OneDrive-YaleUniversity/WCM/OMA/multimedia_data_demo.rda") -tse <- tse_relative +tse <- curatedMetagenomicData( + "HMP_2019_ibdmdb.relative_abundance", + rownames = "short", + dryrun = FALSE +)[[1]] + +# Assign sample IDs for matching. +tse$SampleID <- colnames(tse) + +# Convert relative abundances to the [0, 1] interval. +tse <- transformAssay( + tse, + assay.type = "relative_abundance", + method = "relabundance" +) + +# Use healthy samples as the reference set. +tse$disease_binary <- tse$disease == "healthy" +# Calculate Bray-Curtis dissimilarities across samples. +diss <- as.matrix( + getDissimilarity( + tse, + method = "bray", + na.rm = TRUE, + assay.type = "relabundance" + ) +) +sample_ids <- tse$SampleID + +# Calculate a dysbiosis score for each sample. +tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { + ref_others <- tse$disease_binary & sample_ids != sample_ids[i] + median(diss[i, ref_others], na.rm = TRUE) +}) + +# Keep IBD subjects with both baseline and post-treatment samples. +tse <- tse[, tse$disease != "healthy"] +tse <- tse[, tse$visit_number %in% c(1, 21)] +keep_subjects <- names(which(table(tse$subject_id) == 2)) +tse <- tse[, tse$subject_id %in% keep_subjects] + +# Define treatment as baseline versus post-treatment. +tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") +tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) + +# Standardize species abundances for mediation analysis. +tse <- transformAssay( + tse, + assay.type = "relabundance", + method = "standardize", + name = "scaled" +) +assays(tse) <- assays(tse)["scaled"] ``` The mediation analysis components are defined as following: @@ -208,22 +211,25 @@ treatment's effect on dysbiosis. #| label: overall_mediation_analysis #| message: false -medi_idx <- seq_len(nrow(se_relative)) +# Use indices for mediators. +medi_idx <- seq_len(nrow(tse)) +# Create the mediation data object. exper <- mediation_data( - se_relative, + tse, outcomes = "dysbiosis", treatments = "treatment", mediators = medi_idx ) +# Fit the multimodal mediation model and estimate effects. mdl <- multimedia(exper) res <- estimate(mdl, exper) +# Summarize overall indirect and direct effects. summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) - ``` Mediator-specific indirect effects are then extracted and ranked by absolute @@ -240,7 +246,6 @@ top_mediators <- effects_by_mediator %>% head(20) print(top_mediators) - ``` To quantify uncertainty around the indirect effects, we use non-parametric @@ -256,11 +261,14 @@ boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 10) summary(boot_overall$indirect) quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) -plot_indirect_histogram( - boot_overall$indirect, - "Bootstrap Distribution of Overall Indirect Effect (B = 100)" -) - +ggplot(boot_overall$indirect, aes(indirect_effect)) + + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + + theme_classic() + + labs( + x = "Overall Indirect Effect", + y = "Frequency", + title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" + ) ``` Finally, we visualize the **mediator-specific indirect effects** using a forest @@ -307,11 +315,24 @@ summary_df <- summary_df %>% summary_df$rank <- ifelse(summary_df$significant, 1, 2) summary_df <- summary_df[order(summary_df$rank, -abs(summary_df$estimate)), ] -plot_indirect_forest( - summary_df, - "Forest Plot of Mediation-specific Indirect Effects for Species" -) +summary_df <- summary_df %>% + mutate(abs_estimate = abs(estimate)) +miaViz::plotForest( + summary_df, + effect.var = "estimate", + ci.lower.var = "lower", + ci.upper.var = "upper", + pval.var = NULL, + id.var = "mediator", + order.by = "abs_estimate", + colour.by = "significant" +) + + labs( + x = "Observed Indirect Effect with Bootstrap CI", + y = "Mediator", + title = "Forest Plot of Mediation-specific Indirect Effects for Species" + ) ``` Based on the forest plot results, we can see that the Agathobaculum @@ -339,8 +360,6 @@ ie_pw_fast <- tibble::tibble( ) plot_mediators(ie_pw_fast, exper, n_panels = 12) - - ``` Mediators are standardized (z-scored) for comparability across features; hence values may be negative. Besides, many microbial features are sparse; several mediators exhibit near-zero values for most samples, producing vertical bands. @@ -351,15 +370,62 @@ relative species abundance. ### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} -Let's try to use the pathway abundance data from the iHMP data to do the -analysis as well. +We then load the pathway abundance data, apply the corresponding preprocessing, +and use pathways as mediators. ```{r} #| label: load_iHMP_pathway_demo_data #| message: false -tse <- tse_pathway +tse <- curatedMetagenomicData( + "HMP_2019_ibdmdb.pathway_abundance", + rownames = "short", + dryrun = FALSE +)[[1]] + +# Assign sample IDs for matching. +tse$SampleID <- colnames(tse) + +# Remove samples with missing pathway abundances. +tse <- tse[, colSums(is.na(assay(tse))) == 0] + +# Keep pathway-level features and convert percentages to proportions. +rows_to_keep <- !grepl("\\|", rownames(tse)) +rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE +tse <- tse[rows_to_keep, ] +assay(tse) <- assay(tse) / 100 + +# Use healthy samples as the reference set. +tse$disease_binary <- tse$disease == "healthy" + +# Calculate dysbiosis scores from Bray-Curtis dissimilarities. +diss <- as.matrix(vegdist(t(assay(tse)), method = "bray", na.rm = TRUE)) +tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { + ref_others <- tse$disease_binary & tse$SampleID != tse$SampleID[i] + median(diss[i, ref_others], na.rm = TRUE) +}) + +# Keep IBD subjects with both baseline and post-treatment samples. +tse <- tse[, tse$disease != "healthy"] +tse <- tse[, tse$visit_number %in% c(1, 27)] +keep_subjects <- names(which(table(tse$subject_id) == 2)) +tse <- tse[, tse$subject_id %in% keep_subjects] + +# Define treatment as baseline versus post-treatment. +tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") +tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) + +# Standardize pathway abundances for mediation analysis. +tse <- transformAssay( + tse, + assay.type = "pathway_abundance", + method = "standardize", + name = "scaled" +) +assays(tse) <- assays(tse)["scaled"] +# Clean pathway row names for downstream model terms. +rownames(tse) <- make.names(clean_pathway_names(rownames(tse)), unique = TRUE) ``` Next, we continue to fit the multimodal mediation model, and extract the overall @@ -369,8 +435,10 @@ and mediation-specific indirect effects. #| label: mediation_analysis #| message: false +# Use indices for mediators. medi_idx <- seq_len(nrow(tse)) +# Create the mediation data object. exper <- mediation_data( tse, outcomes = "dysbiosis", @@ -378,8 +446,11 @@ exper <- mediation_data( mediators = medi_idx ) +# Fit the multimodal mediation model and estimate effects. mdl <- multimedia(exper) res <- estimate(mdl, exper) + +# Summarize overall indirect and direct effects. summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) @@ -391,7 +462,6 @@ top_mediators <- effects_by_mediator %>% head(20) print(top_mediators) - ``` The non-parametric bootstrap resampling is also used to compute **95% CIs** for @@ -406,11 +476,14 @@ boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) summary(boot_overall$indirect) quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) -plot_indirect_histogram( - boot_overall$indirect, - "Bootstrap Distribution of Overall Indirect Effect (B = 100)" -) - +ggplot(boot_overall$indirect, aes(indirect_effect)) + + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + + theme_classic() + + labs( + x = "Overall Indirect Effect", + y = "Frequency", + title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" + ) ``` We visualize the **mediator-specific indirect effects** using a forest plot for @@ -455,11 +528,24 @@ pwy <- clean_pathway_names(rownames(summary_df)) rownames(summary_df) <- pwy summary_df$mediator <- pwy -plot_indirect_forest( - summary_df, - "Forest Plot of Mediation-specific Indirect Effects for Pathway" -) +summary_df <- summary_df %>% + mutate(abs_estimate = abs(estimate)) +miaViz::plotForest( + summary_df, + effect.var = "estimate", + ci.lower.var = "lower", + ci.upper.var = "upper", + pval.var = NULL, + id.var = "mediator", + order.by = "abs_estimate", + colour.by = "significant" +) + + labs( + x = "Observed Indirect Effect with Bootstrap CI", + y = "Mediator", + title = "Forest Plot of Mediation-specific Indirect Effects for Pathway" + ) ``` Based on the forest plot results, we can see that the tRNA charging pathway has From 1c45f5c33e7dfbd9bd0160b686f1c094fa5495a5 Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Wed, 22 Jul 2026 02:44:11 +0800 Subject: [PATCH 10/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 238 ++++++++++++++------------------------ 1 file changed, 84 insertions(+), 154 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index f6bad9ab..4739a1b4 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -1,4 +1,4 @@ -# **Multimodal Mediation Analysis** {#sec-MSEA} +# Multimodal Mediation Analysis {#sec-MSEA} ```{r} #| label: setup @@ -30,29 +30,29 @@ layers, providing mechanistic insights into **how exposures impact outcomes** library(DiagrammeR) grViz(" -digraph multimodal_mediation { - graph [layout = dot, rankdir = LR] - - node [fontname = Helvetica, fontsize = 12] - - Exposure [shape=box, style=filled, fillcolor=lightblue, color=black] - Outcome [shape=box, style=filled, fillcolor=palegreen, color=black] - - M1 [label='Species 1', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] - M2 [label='Species 2', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] - M3 [label='Metabolite 1', shape=ellipse, style=filled, fillcolor=lightpink, color=black] - M4 [label='Metabolite 2', shape=ellipse, style=filled, fillcolor=lightpink, color=black] - - Exposure -> M1 - Exposure -> M2 - Exposure -> M3 - Exposure -> M4 - M1 -> Outcome - M2 -> Outcome - M3 -> Outcome - M4 -> Outcome - Exposure -> Outcome [style=dashed] -} + digraph multimodal_mediation { + graph [layout = dot, rankdir = LR] + + node [fontname = Helvetica, fontsize = 12] + + Exposure [shape=box, style=filled, fillcolor=lightblue, color=black] + Outcome [shape=box, style=filled, fillcolor=palegreen, color=black] + + M1 [label='Species 1', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] + M2 [label='Species 2', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] + M3 [label='Metabolite 1', shape=ellipse, style=filled, fillcolor=lightpink, color=black] + M4 [label='Metabolite 2', shape=ellipse, style=filled, fillcolor=lightpink, color=black] + + Exposure -> M1 + Exposure -> M2 + Exposure -> M3 + Exposure -> M4 + M1 -> Outcome + M2 -> Outcome + M3 -> Outcome + M4 -> Outcome + Exposure -> Outcome [style=dashed] + } ") ``` @@ -115,88 +115,30 @@ library(stringr) library(ggplot2) library(mia) library(miaViz) - -clean_pathway_names <- function(pathways) { - pathways <- gsub("PWY0\\.", "PWY0-", pathways) - pathways <- gsub("PWY\\.", "PWY-", pathways) - pathways <- gsub("__", ": ", pathways) - pathways <- gsub("_\\.", " ", pathways) - pathways <- gsub("\\.", " ", pathways) - pathways <- gsub("_", " ", pathways) - pathways <- gsub("\\s+", " ", pathways) - trimws(pathways) -} ``` -We first load the relative abundance data, calculate dysbiosis scores from the -healthy reference samples, and retain IBD subjects with both baseline and -post-treatment samples. +The processed demo data are loaded directly from the local data file prepared +for this chapter. ```{r} #| label: load_iHMP_relative_demo_data #| message: false -tse <- curatedMetagenomicData( - "HMP_2019_ibdmdb.relative_abundance", - rownames = "short", - dryrun = FALSE -)[[1]] +load("multimedia_data_demo.rda") +tse <- tse_relative -# Assign sample IDs for matching. -tse$SampleID <- colnames(tse) - -# Convert relative abundances to the [0, 1] interval. -tse <- transformAssay( - tse, - assay.type = "relative_abundance", - method = "relabundance" -) - -# Use healthy samples as the reference set. -tse$disease_binary <- tse$disease == "healthy" - -# Calculate Bray-Curtis dissimilarities across samples. -diss <- as.matrix( - getDissimilarity( - tse, - method = "bray", - na.rm = TRUE, - assay.type = "relabundance" - ) -) -sample_ids <- tse$SampleID - -# Calculate a dysbiosis score for each sample. -tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { - ref_others <- tse$disease_binary & sample_ids != sample_ids[i] - median(diss[i, ref_others], na.rm = TRUE) -}) - -# Keep IBD subjects with both baseline and post-treatment samples. -tse <- tse[, tse$disease != "healthy"] -tse <- tse[, tse$visit_number %in% c(1, 21)] -keep_subjects <- names(which(table(tse$subject_id) == 2)) -tse <- tse[, tse$subject_id %in% keep_subjects] - -# Define treatment as baseline versus post-treatment. -tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") -tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) - -# Standardize species abundances for mediation analysis. -tse <- transformAssay( - tse, - assay.type = "relabundance", - method = "standardize", - name = "scaled" -) -assays(tse) <- assays(tse)["scaled"] +# Use syntactically valid mediator names in model formulas +raw <- sub(".*s__", "", rownames(tse)) +clean1 <- str_remove_all(raw, "\\[|\\]") +clean2 <- str_replace_all(clean1, "[: \\.,]", "_") +rownames(tse) <- make.names(clean2, unique = TRUE) ``` The mediation analysis components are defined as following: -- Treatment: time point (baseline vs. post-baseline) -- Outcome: dysbiosis score -- Mediators: scaled species abundances +Treatment: time point (baseline vs. post-baseline) +Outcome: dysbiosis score +Mediators: scaled species abundances These are bundled into a mediation data object (exper), which is passed into the multimedia function to estimate path models for the treatment-to-mediator (α @@ -211,10 +153,10 @@ treatment's effect on dysbiosis. #| label: overall_mediation_analysis #| message: false -# Use indices for mediators. +# Use indices for mediators medi_idx <- seq_len(nrow(tse)) -# Create the mediation data object. +# Create the mediation data object exper <- mediation_data( tse, outcomes = "dysbiosis", @@ -222,11 +164,11 @@ exper <- mediation_data( mediators = medi_idx ) -# Fit the multimodal mediation model and estimate effects. +# Fit the multimodal mediation model and estimate effects mdl <- multimedia(exper) res <- estimate(mdl, exper) -# Summarize overall indirect and direct effects. +# Summarize overall indirect and direct effects summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) @@ -241,6 +183,7 @@ effect size. effects_by_mediator <- indirect_pathwise(res, exper) +# Rank mediators by absolute indirect effect size top_mediators <- effects_by_mediator %>% arrange(desc(abs(indirect_effect))) %>% head(20) @@ -256,11 +199,15 @@ overall and mediator-specific indirect effects. #| label: visualization_for_species_overall_indirect_effects #| message: false +# Bootstrap the overall indirect effect set.seed(12345) boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 10) + +# Summarize the bootstrap distribution and percentile interval summary(boot_overall$indirect) quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) +# Plot the bootstrap distribution ggplot(boot_overall$indirect, aes(indirect_effect)) + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + theme_classic() + @@ -279,6 +226,7 @@ mediators. #| label: visualization_for_species_specific_indirect_effects #| message: false +# Define mediator-specific indirect effects for bootstrap resampling indirect_each <- function(mdl, exper) { res <- estimate(mdl, exper) alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) @@ -290,10 +238,12 @@ indirect_each <- function(mdl, exper) { return(indirect) } +# Bootstrap mediator-specific indirect effects set.seed(12345) boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) summary(boot_each$indirect) +# Calculate percentile intervals and bootstrap means lower_upper <- apply( boot_each$indirect, 2, @@ -301,6 +251,7 @@ lower_upper <- apply( ) means <- apply(boot_each$indirect, 2, mean, na.rm = TRUE) +# Collect estimates and interval limits in one data frame summary_df <- data.frame( mediator = colnames(boot_each$indirect), estimate = means, @@ -309,15 +260,22 @@ summary_df <- data.frame( ) summary_df <- summary_df[-1, ] +# Flag intervals that do not cross zero summary_df$significant <- with(summary_df, lower > 0 | upper < 0) + +# Remove mediators with degenerate bootstrap intervals summary_df <- summary_df %>% filter((upper - lower) != 0) + +# Prioritize significant mediators, then sort by effect size summary_df$rank <- ifelse(summary_df$significant, 1, 2) summary_df <- summary_df[order(summary_df$rank, -abs(summary_df$estimate)), ] +# Add absolute effects for forest plot ordering summary_df <- summary_df %>% mutate(abs_estimate = abs(estimate)) +# Plot mediator-specific effects with confidence intervals miaViz::plotForest( summary_df, effect.var = "estimate", @@ -346,11 +304,13 @@ visualization. #| label: default_visualization_plot_mediators_species #| message: false +# Select the top mediators for the default mediator-level plot top_meds <- summary_df %>% dplyr::arrange(dplyr::desc(abs(estimate))) %>% dplyr::slice_head(n = 12) %>% dplyr::pull(mediator) +# Create the effect table expected by plot_mediators ie_pw_fast <- tibble::tibble( outcome = "dysbiosis", mediator = top_meds, @@ -359,6 +319,7 @@ ie_pw_fast <- tibble::tibble( indirect_effect = summary_df$estimate[match(top_meds, summary_df$mediator)] ) +# Visualize the selected mediators across samples plot_mediators(ie_pw_fast, exper, n_panels = 12) ``` Mediators are standardized (z-scored) for comparability across features; hence values may be negative. Besides, many microbial features are sparse; several mediators exhibit near-zero values for most samples, producing vertical bands. @@ -370,62 +331,13 @@ relative species abundance. ### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} -We then load the pathway abundance data, apply the corresponding preprocessing, -and use pathways as mediators. +We then load the processed pathway abundance data and use pathways as mediators. ```{r} #| label: load_iHMP_pathway_demo_data #| message: false -tse <- curatedMetagenomicData( - "HMP_2019_ibdmdb.pathway_abundance", - rownames = "short", - dryrun = FALSE -)[[1]] - -# Assign sample IDs for matching. -tse$SampleID <- colnames(tse) - -# Remove samples with missing pathway abundances. -tse <- tse[, colSums(is.na(assay(tse))) == 0] - -# Keep pathway-level features and convert percentages to proportions. -rows_to_keep <- !grepl("\\|", rownames(tse)) -rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE -tse <- tse[rows_to_keep, ] -assay(tse) <- assay(tse) / 100 - -# Use healthy samples as the reference set. -tse$disease_binary <- tse$disease == "healthy" - -# Calculate dysbiosis scores from Bray-Curtis dissimilarities. -diss <- as.matrix(vegdist(t(assay(tse)), method = "bray", na.rm = TRUE)) -tse$dysbiosis <- sapply(seq_along(tse$disease_binary), function(i) { - ref_others <- tse$disease_binary & tse$SampleID != tse$SampleID[i] - median(diss[i, ref_others], na.rm = TRUE) -}) - -# Keep IBD subjects with both baseline and post-treatment samples. -tse <- tse[, tse$disease != "healthy"] -tse <- tse[, tse$visit_number %in% c(1, 27)] -keep_subjects <- names(which(table(tse$subject_id) == 2)) -tse <- tse[, tse$subject_id %in% keep_subjects] - -# Define treatment as baseline versus post-treatment. -tse$Time_point <- ifelse(tse$visit_number == 1, "T0", "T1") -tse$treatment <- factor(tse$Time_point, levels = c("T0", "T1")) - -# Standardize pathway abundances for mediation analysis. -tse <- transformAssay( - tse, - assay.type = "pathway_abundance", - method = "standardize", - name = "scaled" -) -assays(tse) <- assays(tse)["scaled"] - -# Clean pathway row names for downstream model terms. -rownames(tse) <- make.names(clean_pathway_names(rownames(tse)), unique = TRUE) +tse <- tse_pathway ``` Next, we continue to fit the multimodal mediation model, and extract the overall @@ -435,10 +347,10 @@ and mediation-specific indirect effects. #| label: mediation_analysis #| message: false -# Use indices for mediators. +# Use indices for mediators medi_idx <- seq_len(nrow(tse)) -# Create the mediation data object. +# Create the mediation data object exper <- mediation_data( tse, outcomes = "dysbiosis", @@ -446,17 +358,18 @@ exper <- mediation_data( mediators = medi_idx ) -# Fit the multimodal mediation model and estimate effects. +# Fit the multimodal mediation model and estimate effects mdl <- multimedia(exper) res <- estimate(mdl, exper) -# Summarize overall indirect and direct effects. +# Summarize overall indirect and direct effects summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) effects_by_mediator <- indirect_pathwise(res, exper) +# Rank mediators by absolute indirect effect size top_mediators <- effects_by_mediator %>% arrange(desc(abs(indirect_effect))) %>% head(20) @@ -471,11 +384,15 @@ both the overall and mediator-specific pathway indirect effects. #| label: visualization_for_pathway_overall_indirect_effects #| message: false +# Bootstrap the overall indirect effect set.seed(1234) boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) + +# Summarize the bootstrap distribution and percentile interval summary(boot_overall$indirect) quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) +# Plot the bootstrap distribution ggplot(boot_overall$indirect, aes(indirect_effect)) + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + theme_classic() + @@ -493,6 +410,7 @@ pathway abundance as well, with the point estimates and 95% CIs. #| label: visualization_for_pathway_specific_indirect_effects #| message: false +# Define mediator-specific indirect effects for bootstrap resampling indirect_each <- function(mdl, exper) { res <- estimate(mdl, exper) alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) @@ -504,15 +422,19 @@ indirect_each <- function(mdl, exper) { return(indirect) } +# Bootstrap mediator-specific indirect effects set.seed(1234) boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) +# Calculate percentile intervals and bootstrap means lower_upper <- apply( boot_each$indirect, 2, function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE) ) means <- apply(boot_each$indirect, 2, mean, na.rm = TRUE) + +# Collect estimates and interval limits in one data frame summary_df <- data.frame( mediator = colnames(boot_each$indirect), estimate = means, @@ -521,16 +443,24 @@ summary_df <- data.frame( ) summary_df <- summary_df[-1, ] + +# Flag intervals that do not cross zero summary_df$significant <- with(summary_df, lower > 0 | upper < 0) + +# Remove mediators with degenerate bootstrap intervals summary_df <- summary_df[summary_df$upper != summary_df$lower, ] -pwy <- clean_pathway_names(rownames(summary_df)) +# Restore readable pathway labels for plotting +pwy <- gsub("PWY0\\.|PWY\\.", "PWY-", rownames(summary_df)) +pwy <- trimws(gsub("\\s+", " ", gsub("[_\\.]+", " ", pwy))) rownames(summary_df) <- pwy summary_df$mediator <- pwy +# Add absolute effects for forest plot ordering summary_df <- summary_df %>% mutate(abs_estimate = abs(estimate)) +# Plot mediator-specific effects with confidence intervals miaViz::plotForest( summary_df, effect.var = "estimate", From 1d7d75c25cbb579877aa915289a3f7446b37ca4f Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:35:30 +0800 Subject: [PATCH 11/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 189 +++++++++++++++++++++++++++----------- 1 file changed, 137 insertions(+), 52 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 4739a1b4..1bb6bc14 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -27,33 +27,44 @@ layers, providing mechanistic insights into **how exposures impact outcomes** #| fig-cap: Directed acyclic graph illustrating multimodal mediation where an exposure affects multiple mediators across modalities (microbial species and metabolites), which in turn affect the outcome. A direct path from exposure to outcome is also included (dashed). While gut microbial species biologically influence metabolite levels, for tractability in high-dimensional mediation analysis we assume species and metabolites act as conditionally independent parallel mediators given exposure, acknowledging this as an approximation #| echo: false -library(DiagrammeR) - -grViz(" - digraph multimodal_mediation { - graph [layout = dot, rankdir = LR] - - node [fontname = Helvetica, fontsize = 12] - - Exposure [shape=box, style=filled, fillcolor=lightblue, color=black] - Outcome [shape=box, style=filled, fillcolor=palegreen, color=black] - - M1 [label='Species 1', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] - M2 [label='Species 2', shape=ellipse, style=filled, fillcolor=lightyellow, color=black] - M3 [label='Metabolite 1', shape=ellipse, style=filled, fillcolor=lightpink, color=black] - M4 [label='Metabolite 2', shape=ellipse, style=filled, fillcolor=lightpink, color=black] - - Exposure -> M1 - Exposure -> M2 - Exposure -> M3 - Exposure -> M4 - M1 -> Outcome - M2 -> Outcome - M3 -> Outcome - M4 -> Outcome - Exposure -> Outcome [style=dashed] - } -") +library(ggdag) +library(ggplot2) + +# Define DAG: exposure -> 4 parallel mediators -> outcome, plus a direct path +dag <- dagify( + outcome ~ sp1 + sp2 + met1 + met2 + exposure, + sp1 ~ exposure, + sp2 ~ exposure, + met1 ~ exposure, + met2 ~ exposure, + exposure = "exposure", + outcome = "outcome", + labels = c( + exposure = "Exposure", + outcome = "Outcome", + sp1 = "Species 1", + sp2 = "Species 2", + met1 = "Metabolite 1", + met2 = "Metabolite 2" + ), + coords = list( + x = c(exposure = 0, sp1 = 1, sp2 = 1, met1 = 1, met2 = 1, outcome = 2), + y = c(exposure = 0, sp1 = 1.5, sp2 = 0.5, met1 = -0.5, met2 = -1.5, outcome = 0) + ) +) + +# Flag the direct exposure -> outcome edge so it renders dashed +tidy_dag <- tidy_dagitty(dag) +tidy_dag$data$direct <- with( + tidy_dag$data, name == "exposure" & !is.na(to) & to == "outcome" +) + +ggplot(tidy_dag, aes(x = x, y = y, xend = xend, yend = yend)) + + geom_dag_edges(aes(edge_linetype = direct)) + + geom_dag_point(color = "lightblue") + + geom_dag_text(aes(label = label), color = "black", size = 3) + + scale_edge_linetype_manual(values = c("solid", "dashed"), guide = "none") + + theme_dag() ``` In this chapter, we demonstrate multimodal mediation analysis using microbiome @@ -91,7 +102,7 @@ prioritize findings. 7. Repeat the process for the iHMP microbial pathways. -### Performing multimodal mediation analysis for iHMP species relative abundance {#sec-relative-abundance} +## Performing multimodal mediation analysis for iHMP species relative abundance {#sec-relative-abundance} We begin by loading the relative abundance data and calculating dysbiosis scores for each sample based on Bray–Curtis dissimilarity from a healthy reference set. @@ -117,21 +128,56 @@ library(mia) library(miaViz) ``` -The processed demo data are loaded directly from the local data file prepared -for this chapter. +We pull the iHMP IBDMDB species relative-abundance table directly from +`curatedMetagenomicData` and apply a minimal preprocessing pipeline: convert +to relative abundances, compute per-sample dysbiosis scores (median +Bray–Curtis distance from the healthy reference), subset to IBD subjects +with paired baseline (visit 1) and post-treatment (visit 21) samples, +define the treatment factor, and standardize the assay for mediation. ```{r} #| label: load_iHMP_relative_demo_data #| message: false -load("multimedia_data_demo.rda") -tse <- tse_relative +# Import iHMP species relative abundance data +tse <- curatedMetagenomicData( + "HMP_2019_ibdmdb.relative_abundance", + rownames = "short", dryrun = FALSE +)[[1]] + +# Convert the assay to relative abundances in [0, 1] +tse <- transformAssay( + tse, assay.type = "relative_abundance", method = "relabundance" +) + +# Dysbiosis: median Bray-Curtis distance to healthy references (excluding self) +diss <- as.matrix(getDissimilarity( + tse, method = "bray", assay.type = "relabundance", na.rm = TRUE +)) +healthy <- tse$disease == "healthy" +tse$dysbiosis <- sapply(seq_len(ncol(tse)), function(i) { + median(diss[i, healthy & seq_len(ncol(tse)) != i], na.rm = TRUE) +}) + +# Keep IBD subjects with paired baseline (visit 1) and post-treatment (visit 21) +tse <- tse[, tse$disease != "healthy" & tse$visit_number %in% c(1, 21)] +paired <- names(which(table(tse$subject_id) == 2)) +tse <- tse[, tse$subject_id %in% paired] + +# Treatment factor and standardized mediator assay +tse$treatment <- factor( + ifelse(tse$visit_number == 1, "T0", "T1"), levels = c("T0", "T1") +) +tse <- transformAssay( + tse, assay.type = "relabundance", method = "standardize", name = "scaled" +) # Use syntactically valid mediator names in model formulas -raw <- sub(".*s__", "", rownames(tse)) -clean1 <- str_remove_all(raw, "\\[|\\]") -clean2 <- str_replace_all(clean1, "[: \\.,]", "_") -rownames(tse) <- make.names(clean2, unique = TRUE) +rownames(tse) <- rownames(tse) |> + str_remove(".*s__") |> + str_remove_all("\\[|\\]") |> + str_replace_all("[: \\.,]", "_") |> + make.names(unique = TRUE) ``` The mediation analysis components are defined as following: @@ -276,7 +322,7 @@ summary_df <- summary_df %>% mutate(abs_estimate = abs(estimate)) # Plot mediator-specific effects with confidence intervals -miaViz::plotForest( +plotForest( summary_df, effect.var = "estimate", ci.lower.var = "lower", @@ -306,12 +352,12 @@ visualization. # Select the top mediators for the default mediator-level plot top_meds <- summary_df %>% - dplyr::arrange(dplyr::desc(abs(estimate))) %>% - dplyr::slice_head(n = 12) %>% - dplyr::pull(mediator) + arrange(desc(abs(estimate))) %>% + slice_head(n = 12) %>% + pull(mediator) # Create the effect table expected by plot_mediators -ie_pw_fast <- tibble::tibble( +ie_pw_fast <- tibble( outcome = "dysbiosis", mediator = top_meds, direct_setting = levels(tse$treatment)[1], @@ -329,15 +375,54 @@ same setup as before, using the time point as the treatment, dysbiosis score as the outcome, but pathways abundances as the mediation variables instead of the relative species abundance. -### Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} +## Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} -We then load the processed pathway abundance data and use pathways as mediators. +We then repeat the same minimal preprocessing on the iHMP pathway abundance +table (paired visits 1 and 27 for pathways). ```{r} #| label: load_iHMP_pathway_demo_data #| message: false -tse <- tse_pathway +# Import iHMP pathway abundance data +tse <- curatedMetagenomicData( + "HMP_2019_ibdmdb.pathway_abundance", + rownames = "short", dryrun = FALSE +)[[1]] + +# Drop samples with any NAs; keep only top-level pathways (rows without "|" +# strata) and remove the UNMAPPED / UNINTEGRATED rows; scale counts to [0, 1] +tse <- tse[, colSums(is.na(assay(tse))) == 0] +rows_to_keep <- !grepl("\\|", rownames(tse)) +rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE +tse <- tse[rows_to_keep, ] +assay(tse) <- assay(tse) / 100 + +# Dysbiosis: median Bray-Curtis distance to healthy references (excluding self) +diss <- as.matrix(vegdist(t(assay(tse)), method = "bray", na.rm = TRUE)) +healthy <- tse$disease == "healthy" +tse$dysbiosis <- sapply(seq_len(ncol(tse)), function(i) { + median(diss[i, healthy & seq_len(ncol(tse)) != i], na.rm = TRUE) +}) + +# Keep IBD subjects with paired baseline (visit 1) and post-treatment (visit 27) +tse <- tse[, tse$disease != "healthy" & tse$visit_number %in% c(1, 27)] +paired <- names(which(table(tse$subject_id) == 2)) +tse <- tse[, tse$subject_id %in% paired] + +# Treatment factor and standardized mediator assay +tse$treatment <- factor( + ifelse(tse$visit_number == 1, "T0", "T1"), levels = c("T0", "T1") +) +tse <- transformAssay( + tse, assay.type = "pathway_abundance", method = "standardize", name = "scaled" +) + +# Syntactically valid mediator names for model formulas +rownames(tse) <- rownames(tse) |> + str_remove_all("\\[|\\]") |> + str_replace_all("[: \\.,]", "_") |> + make.names(unique = TRUE) ``` Next, we continue to fit the multimodal mediation model, and extract the overall @@ -435,11 +520,12 @@ lower_upper <- apply( means <- apply(boot_each$indirect, 2, mean, na.rm = TRUE) # Collect estimates and interval limits in one data frame +# (plotForest reads rownames by default, so no mediator column needed) summary_df <- data.frame( - mediator = colnames(boot_each$indirect), estimate = means, lower = lower_upper[1, ], - upper = lower_upper[2, ] + upper = lower_upper[2, ], + row.names = colnames(boot_each$indirect) ) summary_df <- summary_df[-1, ] @@ -451,23 +537,22 @@ summary_df$significant <- with(summary_df, lower > 0 | upper < 0) summary_df <- summary_df[summary_df$upper != summary_df$lower, ] # Restore readable pathway labels for plotting -pwy <- gsub("PWY0\\.|PWY\\.", "PWY-", rownames(summary_df)) -pwy <- trimws(gsub("\\s+", " ", gsub("[_\\.]+", " ", pwy))) -rownames(summary_df) <- pwy -summary_df$mediator <- pwy +rownames(summary_df) <- rownames(summary_df) |> + str_replace_all("PWY0\\.|PWY\\.", "PWY-") |> + str_replace_all("[_\\.]+", " ") |> + str_squish() # Add absolute effects for forest plot ordering summary_df <- summary_df %>% mutate(abs_estimate = abs(estimate)) # Plot mediator-specific effects with confidence intervals -miaViz::plotForest( +plotForest( summary_df, effect.var = "estimate", ci.lower.var = "lower", ci.upper.var = "upper", pval.var = NULL, - id.var = "mediator", order.by = "abs_estimate", colour.by = "significant" ) + From c5df3c07686306d1710591fdb2081be7be8cdf7d Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Wed, 5 Aug 2026 03:15:15 +0800 Subject: [PATCH 12/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 1bb6bc14..fb981c67 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -147,12 +147,14 @@ tse <- curatedMetagenomicData( # Convert the assay to relative abundances in [0, 1] tse <- transformAssay( - tse, assay.type = "relative_abundance", method = "relabundance" + tse, + assay.type = "relative_abundance", method = "relabundance" ) # Dysbiosis: median Bray-Curtis distance to healthy references (excluding self) diss <- as.matrix(getDissimilarity( - tse, method = "bray", assay.type = "relabundance", na.rm = TRUE + tse, + method = "bray", assay.type = "relabundance", na.rm = TRUE )) healthy <- tse$disease == "healthy" tse$dysbiosis <- sapply(seq_len(ncol(tse)), function(i) { @@ -166,10 +168,12 @@ tse <- tse[, tse$subject_id %in% paired] # Treatment factor and standardized mediator assay tse$treatment <- factor( - ifelse(tse$visit_number == 1, "T0", "T1"), levels = c("T0", "T1") + ifelse(tse$visit_number == 1, "T0", "T1"), + levels = c("T0", "T1") ) tse <- transformAssay( - tse, assay.type = "relabundance", method = "standardize", name = "scaled" + tse, + assay.type = "relabundance", method = "standardize", name = "scaled" ) # Use syntactically valid mediator names in model formulas @@ -412,10 +416,12 @@ tse <- tse[, tse$subject_id %in% paired] # Treatment factor and standardized mediator assay tse$treatment <- factor( - ifelse(tse$visit_number == 1, "T0", "T1"), levels = c("T0", "T1") + ifelse(tse$visit_number == 1, "T0", "T1"), + levels = c("T0", "T1") ) tse <- transformAssay( - tse, assay.type = "pathway_abundance", method = "standardize", name = "scaled" + tse, + assay.type = "pathway_abundance", method = "standardize", name = "scaled" ) # Syntactically valid mediator names for model formulas From 9899618449f9c79643b722999c531d8af557fac4 Mon Sep 17 00:00:00 2001 From: Giulio Date: Thu, 6 Aug 2026 16:06:34 +0300 Subject: [PATCH 13/16] Refine multimedia chapter --- inst/pages/multimedia.qmd | 284 ++++++++++++++++++++------------------ 1 file changed, 148 insertions(+), 136 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index fb981c67..be823d44 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -28,6 +28,7 @@ layers, providing mechanistic insights into **how exposures impact outcomes** #| echo: false library(ggdag) +library(ggraph) library(ggplot2) # Define DAG: exposure -> 4 parallel mediators -> outcome, plus a direct path @@ -116,16 +117,14 @@ We first load the packages and define small plotting helpers used below. #| label: load_libraries_and_helpers #| message: false +library(mia) +library(miaViz) +library(multimedia) library(curatedMetagenomicData) -library(SummarizedExperiment) library(dplyr) -library(vegan) library(tidyverse) -library(multimedia) library(stringr) library(ggplot2) -library(mia) -library(miaViz) ``` We pull the iHMP IBDMDB species relative-abundance table directly from @@ -139,30 +138,39 @@ define the treatment factor, and standardize the assay for mediation. #| label: load_iHMP_relative_demo_data #| message: false -# Import iHMP species relative abundance data -tse <- curatedMetagenomicData( +# Import iHMP dataset +hmpibd <- curatedMetagenomicData( "HMP_2019_ibdmdb.relative_abundance", - rownames = "short", dryrun = FALSE -)[[1]] + rownames = "short", + dryrun = FALSE +) + +# Extract microbiome experiment +tse <- hmpibd[[1]] # Convert the assay to relative abundances in [0, 1] tse <- transformAssay( tse, - assay.type = "relative_abundance", method = "relabundance" + assay.type = "relative_abundance", + method = "relabundance" ) # Dysbiosis: median Bray-Curtis distance to healthy references (excluding self) diss <- as.matrix(getDissimilarity( tse, - method = "bray", assay.type = "relabundance", na.rm = TRUE + method = "bray", + assay.type = "relabundance", + na.rm = TRUE )) -healthy <- tse$disease == "healthy" + +is_healthy <- tse$disease == "healthy" + tse$dysbiosis <- sapply(seq_len(ncol(tse)), function(i) { - median(diss[i, healthy & seq_len(ncol(tse)) != i], na.rm = TRUE) + median(diss[i, is_healthy & seq_len(ncol(tse)) != i], na.rm = TRUE) }) # Keep IBD subjects with paired baseline (visit 1) and post-treatment (visit 21) -tse <- tse[, tse$disease != "healthy" & tse$visit_number %in% c(1, 21)] +tse <- tse[, !is_healthy & tse$visit_number %in% c(1, 21)] paired <- names(which(table(tse$subject_id) == 2)) tse <- tse[, tse$subject_id %in% paired] @@ -171,9 +179,12 @@ tse$treatment <- factor( ifelse(tse$visit_number == 1, "T0", "T1"), levels = c("T0", "T1") ) + tse <- transformAssay( tse, - assay.type = "relabundance", method = "standardize", name = "scaled" + assay.type = "relabundance", + method = "standardize", + name = "scaled" ) # Use syntactically valid mediator names in model formulas @@ -234,8 +245,8 @@ effect size. effects_by_mediator <- indirect_pathwise(res, exper) # Rank mediators by absolute indirect effect size -top_mediators <- effects_by_mediator %>% - arrange(desc(abs(indirect_effect))) %>% +top_mediators <- effects_by_mediator |> + arrange(desc(abs(indirect_effect))) |> head(20) print(top_mediators) @@ -249,8 +260,9 @@ overall and mediator-specific indirect effects. #| label: visualization_for_species_overall_indirect_effects #| message: false -# Bootstrap the overall indirect effect set.seed(12345) + +# Bootstrap the overall indirect effect boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 10) # Summarize the bootstrap distribution and percentile interval @@ -260,12 +272,9 @@ quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) # Plot the bootstrap distribution ggplot(boot_overall$indirect, aes(indirect_effect)) + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + - theme_classic() + - labs( - x = "Overall Indirect Effect", - y = "Frequency", - title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" - ) + ggtitle("Bootstrap Distribution of Overall Indirect Effect (B = 100)") + + labs(x = "Overall Indirect Effect", y = "Frequency") + + theme_classic() ``` Finally, we visualize the **mediator-specific indirect effects** using a forest @@ -288,59 +297,53 @@ indirect_each <- function(mdl, exper) { return(indirect) } -# Bootstrap mediator-specific indirect effects set.seed(12345) + +# Bootstrap mediator-specific indirect effects boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) -summary(boot_each$indirect) + +# Remove bootstrap index +boot_each <- boot_each$indirect[, -1] + +summary(boot_each) # Calculate percentile intervals and bootstrap means lower_upper <- apply( - boot_each$indirect, - 2, - function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE) + boot_each, + MARGIN = 2, + quantile, + probs = c(0.025, 0.975), + na.rm = TRUE ) -means <- apply(boot_each$indirect, 2, mean, na.rm = TRUE) # Collect estimates and interval limits in one data frame summary_df <- data.frame( - mediator = colnames(boot_each$indirect), - estimate = means, + mediator = colnames(boot_each), lower = lower_upper[1, ], upper = lower_upper[2, ] ) -summary_df <- summary_df[-1, ] - -# Flag intervals that do not cross zero -summary_df$significant <- with(summary_df, lower > 0 | upper < 0) -# Remove mediators with degenerate bootstrap intervals -summary_df <- summary_df %>% - filter((upper - lower) != 0) +summary_df$effect <- apply(boot_each, 2, mean, na.rm = TRUE) -# Prioritize significant mediators, then sort by effect size -summary_df$rank <- ifelse(summary_df$significant, 1, 2) -summary_df <- summary_df[order(summary_df$rank, -abs(summary_df$estimate)), ] - -# Add absolute effects for forest plot ordering -summary_df <- summary_df %>% - mutate(abs_estimate = abs(estimate)) +summary_df <- summary_df |> + mutate( + # Flag intervals that do not cross zero + significant = lower > 0 | upper < 0, + # Add absolute effects for forest plot ordering + abs_effect = abs(effect) + ) |> + # Remove mediators with degenerate bootstrap intervals + filter(upper - lower > 0) # Plot mediator-specific effects with confidence intervals plotForest( summary_df, - effect.var = "estimate", - ci.lower.var = "lower", - ci.upper.var = "upper", - pval.var = NULL, id.var = "mediator", - order.by = "abs_estimate", - colour.by = "significant" + label.by = "CI", + order.by = "abs_effect" ) + - labs( - x = "Observed Indirect Effect with Bootstrap CI", - y = "Mediator", - title = "Forest Plot of Mediation-specific Indirect Effects for Species" - ) + ggtitle("Forest Plot of Mediation-specific Indirect Effects for Species") + + labs(x = "Observed Indirect Effect with Bootstrap CI", y = "Mediator") ``` Based on the forest plot results, we can see that the Agathobaculum @@ -355,9 +358,9 @@ visualization. #| message: false # Select the top mediators for the default mediator-level plot -top_meds <- summary_df %>% - arrange(desc(abs(estimate))) %>% - slice_head(n = 12) %>% +top_meds <- summary_df |> + arrange(desc(abs_effect)) |> + slice_head(n = 12) |> pull(mediator) # Create the effect table expected by plot_mediators @@ -372,6 +375,7 @@ ie_pw_fast <- tibble( # Visualize the selected mediators across samples plot_mediators(ie_pw_fast, exper, n_panels = 12) ``` + Mediators are standardized (z-scored) for comparability across features; hence values may be negative. Besides, many microbial features are sparse; several mediators exhibit near-zero values for most samples, producing vertical bands. Next we repeat the mediation analysis with the iHMP pathway abundances using the @@ -389,43 +393,63 @@ table (paired visits 1 and 27 for pathways). #| message: false # Import iHMP pathway abundance data -tse <- curatedMetagenomicData( +hmpibd <- curatedMetagenomicData( "HMP_2019_ibdmdb.pathway_abundance", - rownames = "short", dryrun = FALSE -)[[1]] + rownames = "short", + dryrun = FALSE +) + +# Extract pathway experiment from list +se <- hmpibd[[1]] + +# Drop samples with any NAs +se <- se[, colSums(is.na(assay(se))) == 0] -# Drop samples with any NAs; keep only top-level pathways (rows without "|" -# strata) and remove the UNMAPPED / UNINTEGRATED rows; scale counts to [0, 1] -tse <- tse[, colSums(is.na(assay(tse))) == 0] -rows_to_keep <- !grepl("\\|", rownames(tse)) -rows_to_keep[which(rows_to_keep)[1:2]] <- FALSE -tse <- tse[rows_to_keep, ] -assay(tse) <- assay(tse) / 100 +# Keep only top-level pathways (rows without "|" strata) +se <- se[str_detect(rownames(se), fixed("|"), negate = TRUE), ] + +# Remove the UNMAPPED / UNINTEGRATED rows +se <- se[str_detect(rownames(se), "UNMAPPED|UNINTEGRATED", negate = TRUE), ] + +# Scale counts to [0, 1] +assay(se) <- assay(se) / 100 # Dysbiosis: median Bray-Curtis distance to healthy references (excluding self) -diss <- as.matrix(vegdist(t(assay(tse)), method = "bray", na.rm = TRUE)) -healthy <- tse$disease == "healthy" -tse$dysbiosis <- sapply(seq_len(ncol(tse)), function(i) { - median(diss[i, healthy & seq_len(ncol(tse)) != i], na.rm = TRUE) +diss <- as.matrix(getDissimilarity( + se, + method = "bray", + assay.type = "pathway_abundance", + na.rm = TRUE +)) + +is_healthy <- se$disease == "healthy" + +se$dysbiosis <- sapply(seq_len(ncol(se)), function(i) { + median(diss[i, is_healthy & seq_len(ncol(se)) != i], na.rm = TRUE) }) # Keep IBD subjects with paired baseline (visit 1) and post-treatment (visit 27) -tse <- tse[, tse$disease != "healthy" & tse$visit_number %in% c(1, 27)] -paired <- names(which(table(tse$subject_id) == 2)) -tse <- tse[, tse$subject_id %in% paired] +se <- se[, !is_healthy & se$visit_number %in% c(1, 27)] + +paired <- names(which(table(se$subject_id) == 2)) + +se <- se[, se$subject_id %in% paired] # Treatment factor and standardized mediator assay -tse$treatment <- factor( - ifelse(tse$visit_number == 1, "T0", "T1"), +se$treatment <- factor( + ifelse(se$visit_number == 1, "T0", "T1"), levels = c("T0", "T1") ) -tse <- transformAssay( - tse, - assay.type = "pathway_abundance", method = "standardize", name = "scaled" + +se <- transformAssay( + se, + assay.type = "pathway_abundance", + method = "standardize", + name = "scaled" ) # Syntactically valid mediator names for model formulas -rownames(tse) <- rownames(tse) |> +rownames(se) <- rownames(se) |> str_remove_all("\\[|\\]") |> str_replace_all("[: \\.,]", "_") |> make.names(unique = TRUE) @@ -439,11 +463,11 @@ and mediation-specific indirect effects. #| message: false # Use indices for mediators -medi_idx <- seq_len(nrow(tse)) +medi_idx <- seq_len(nrow(se)) # Create the mediation data object exper <- mediation_data( - tse, + se, outcomes = "dysbiosis", treatments = "treatment", mediators = medi_idx @@ -457,12 +481,15 @@ res <- estimate(mdl, exper) summary(res) print(indirect_overall(res, exper)) print(direct_effect(res, exper)) +``` + +```{r} effects_by_mediator <- indirect_pathwise(res, exper) # Rank mediators by absolute indirect effect size -top_mediators <- effects_by_mediator %>% - arrange(desc(abs(indirect_effect))) %>% +top_mediators <- effects_by_mediator |> + arrange(desc(abs(indirect_effect))) |> head(20) print(top_mediators) @@ -475,8 +502,9 @@ both the overall and mediator-specific pathway indirect effects. #| label: visualization_for_pathway_overall_indirect_effects #| message: false +set.seed(12345) + # Bootstrap the overall indirect effect -set.seed(1234) boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) # Summarize the bootstrap distribution and percentile interval @@ -486,12 +514,9 @@ quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) # Plot the bootstrap distribution ggplot(boot_overall$indirect, aes(indirect_effect)) + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + - theme_classic() + - labs( - x = "Overall Indirect Effect", - y = "Frequency", - title = "Bootstrap Distribution of Overall Indirect Effect (B = 100)" - ) + ggtitle("Bootstrap Distribution of Overall Indirect Effect (B = 100)") + + labs(x = "Overall Indirect Effect", y = "Frequency") + + theme_classic() ``` We visualize the **mediator-specific indirect effects** using a forest plot for @@ -501,72 +526,59 @@ pathway abundance as well, with the point estimates and 95% CIs. #| label: visualization_for_pathway_specific_indirect_effects #| message: false -# Define mediator-specific indirect effects for bootstrap resampling -indirect_each <- function(mdl, exper) { - res <- estimate(mdl, exper) - alpha <- sapply(res@mediation@estimates, function(m) coef(m)["treatmentT1"]) - names(alpha) <- sub("\\.treatmentT1$", "", names(alpha)) - beta <- coef(res@outcome@estimates$dysbiosis)[names(alpha)] - indirect <- alpha * beta - names(indirect) <- names(alpha) - indirect <- indirect[!is.na(indirect)] - return(indirect) -} +set.seed(12345) # Bootstrap mediator-specific indirect effects -set.seed(1234) boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) +# Remove bootstrap index +boot_each <- boot_each$indirect[, -1] + +summary(boot_each) + # Calculate percentile intervals and bootstrap means lower_upper <- apply( - boot_each$indirect, - 2, - function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE) + boot_each, + MARGIN = 2, + quantile, + probs = c(0.025, 0.975), + na.rm = TRUE ) -means <- apply(boot_each$indirect, 2, mean, na.rm = TRUE) # Collect estimates and interval limits in one data frame -# (plotForest reads rownames by default, so no mediator column needed) summary_df <- data.frame( - estimate = means, + mediator = colnames(boot_each), lower = lower_upper[1, ], - upper = lower_upper[2, ], - row.names = colnames(boot_each$indirect) + upper = lower_upper[2, ] ) -summary_df <- summary_df[-1, ] +summary_df$effect <- apply(boot_each, 2, mean, na.rm = TRUE) -# Flag intervals that do not cross zero -summary_df$significant <- with(summary_df, lower > 0 | upper < 0) - -# Remove mediators with degenerate bootstrap intervals -summary_df <- summary_df[summary_df$upper != summary_df$lower, ] +summary_df <- summary_df |> + mutate( + # Flag intervals that do not cross zero + significant = lower > 0 | upper < 0, + # Add absolute effects for forest plot ordering + abs_effect = abs(effect) + ) |> + # Remove mediators with degenerate bootstrap intervals + filter(upper - lower > 0) # Restore readable pathway labels for plotting -rownames(summary_df) <- rownames(summary_df) |> +summary_df$mediator <- summary_df$mediator |> str_replace_all("PWY0\\.|PWY\\.", "PWY-") |> str_replace_all("[_\\.]+", " ") |> str_squish() -# Add absolute effects for forest plot ordering -summary_df <- summary_df %>% - mutate(abs_estimate = abs(estimate)) - # Plot mediator-specific effects with confidence intervals plotForest( summary_df, - effect.var = "estimate", - ci.lower.var = "lower", - ci.upper.var = "upper", - pval.var = NULL, - order.by = "abs_estimate", - colour.by = "significant" + id.var = "mediator", + label.by = "CI", + order.by = "abs_effect" ) + - labs( - x = "Observed Indirect Effect with Bootstrap CI", - y = "Mediator", - title = "Forest Plot of Mediation-specific Indirect Effects for Pathway" - ) + ggtitle("Forest Plot of Mediation-specific Indirect Effects for Pathway") + + labs(x = "Observed Indirect Effect with Bootstrap CI", y = "Mediator") ``` Based on the forest plot results, we can see that the tRNA charging pathway has From c8a3a5533a9c9a89caaf3db000455bf371283156 Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Sun, 9 Aug 2026 13:07:04 +0800 Subject: [PATCH 14/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 95 ++++++++++++++++++++------------------- 1 file changed, 48 insertions(+), 47 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index be823d44..6c2b98e8 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -1,4 +1,4 @@ -# Multimodal Mediation Analysis {#sec-MSEA} +# Multivariate Mediation Analysis {#sec-MSEA} ```{r} #| label: setup @@ -10,21 +10,23 @@ library(rebook) chapterPreamble() ``` -Building upon the concepts and workflows presented in the previous chapter on -mediation analysis, here we demonstrate how to perform **multimodal mediation** -**analysis** to identify potential mediators across multiple omics layers. -**Multimodal mediation analysis** examines whether and how an exposure (e.g., -treatment, diet) affects an outcome (e.g., disease, phenotype) **through** -**intermediate variables (mediators) across multiple omics layers** -**simultaneously** (e.g., microbiome, metabolomics). It identifies which -features across modalities act as mediators, quantifies indirect (mediated) -effects while adjusting for covariates, and compares mediation strength across -layers, providing mechanistic insights into **how exposures impact outcomes** -**via molecular pathways**. +Building upon the concepts and workflows presented in the previous chapter on +mediation analysis, here we demonstrate how to perform **multivariate** +**mediation analysis** to identify potential mediators when many candidate +features are measured simultaneously within a single omic layer. +**Multivariate mediation analysis** examines whether and how an exposure +(e.g., treatment, diet) affects an outcome (e.g., disease, phenotype) +**through a high-dimensional set of intermediate variables (mediators)**, +such as the hundreds of species-level abundances or pathway abundances +in a microbiome study. It identifies which features act as mediators, +quantifies indirect (mediated) effects while adjusting for covariates, +and ranks mediation strength across features, providing mechanistic +insights into **how exposures impact outcomes via microbial features and** +**functional pathways**. ```{r} -#| label: fig_multimodal_mediation -#| fig-cap: Directed acyclic graph illustrating multimodal mediation where an exposure affects multiple mediators across modalities (microbial species and metabolites), which in turn affect the outcome. A direct path from exposure to outcome is also included (dashed). While gut microbial species biologically influence metabolite levels, for tractability in high-dimensional mediation analysis we assume species and metabolites act as conditionally independent parallel mediators given exposure, acknowledging this as an approximation +#| label: fig_multivariate_mediation +#| fig-cap: Directed acyclic graph illustrating multivariate mediation, where an exposure affects many parallel mediators (here, illustrated as microbial features) within a single omic layer, which in turn affect the outcome. A direct path from exposure to outcome is also included (dashed). For tractability in high-dimensional mediation analysis, mediators are assumed to be conditionally independent given the exposure; this is a simplification, since real microbial features are often correlated. #| echo: false library(ggdag) @@ -33,24 +35,24 @@ library(ggplot2) # Define DAG: exposure -> 4 parallel mediators -> outcome, plus a direct path dag <- dagify( - outcome ~ sp1 + sp2 + met1 + met2 + exposure, - sp1 ~ exposure, - sp2 ~ exposure, - met1 ~ exposure, - met2 ~ exposure, + outcome ~ m1 + m2 + m3 + m4 + exposure, + m1 ~ exposure, + m2 ~ exposure, + m3 ~ exposure, + m4 ~ exposure, exposure = "exposure", outcome = "outcome", labels = c( exposure = "Exposure", outcome = "Outcome", - sp1 = "Species 1", - sp2 = "Species 2", - met1 = "Metabolite 1", - met2 = "Metabolite 2" + m1 = "Feature 1", + m2 = "Feature 2", + m3 = "Feature 3", + m4 = "Feature 4" ), coords = list( - x = c(exposure = 0, sp1 = 1, sp2 = 1, met1 = 1, met2 = 1, outcome = 2), - y = c(exposure = 0, sp1 = 1.5, sp2 = 0.5, met1 = -0.5, met2 = -1.5, outcome = 0) + x = c(exposure = 0, m1 = 1, m2 = 1, m3 = 1, m4 = 1, outcome = 2), + y = c(exposure = 0, m1 = 1.5, m2 = 0.5, m3 = -0.5, m4 = -1.5, outcome = 0) ) ) @@ -68,17 +70,16 @@ ggplot(tidy_dag, aes(x = x, y = y, xend = xend, yend = yend)) + theme_dag() ``` -In this chapter, we demonstrate multimodal mediation analysis using microbiome -data from the iHMP IBDMDB study from the R/Bioconductor package -curatedMetagenomicData. [@Lloyd-Price2019]. Unlike conventional mediation -analysis, which typically focuses on a single mediator or a small set of -mediators, we performed multimodal mediation analysis with R/Bioconductor -package \[multimedia\] [@Jiang2025] here, which can handle many potential -mediators across multiple data modalities. For example, species-level taxonomic -abundances, pathways abundances, and metabolomic profiles can all serve as -potential mediators. This method allows us to estimate and rank each mediator’s -contribution while accounting for correlations between mediators in a high- -dimensional setting. +In this chapter, we demonstrate multivariate mediation analysis using microbiome +data from the iHMP IBDMDB study from the R/Bioconductor package +curatedMetagenomicData. [@Lloyd-Price2019]. Unlike conventional mediation +analysis, which typically focuses on a single mediator or a small set of +mediators, we performed multivariate mediation analysis with the R/Bioconductor +package \[multimedia\] [@Jiang2025] here, which can handle many potential +mediators within a single omic layer. For example, species-level taxonomic +abundances or pathway abundances can each serve as a set of parallel mediators. +This method allows us to estimate and rank each mediator’s contribution while +accounting for correlations between mediators in a high-dimensional setting. In this example, time point (baseline vs. post-treatment) is used as the treatment, and microbiome-based dysbiosis scores are used as the outcome. We @@ -93,7 +94,7 @@ Generally, we will proceed through the following key steps: 3. Define the mediation analysis data structure. -4. Fit the multimodal mediation model using the R package multimedia. +4. Fit the multivariate mediation model using the R package multimedia. 5. Interpret both the overall indirect effects and the mediator-specific indirect effects. @@ -103,7 +104,7 @@ prioritize findings. 7. Repeat the process for the iHMP microbial pathways. -## Performing multimodal mediation analysis for iHMP species relative abundance {#sec-relative-abundance} +## Performing multivariate mediation analysis for iHMP species relative abundance {#sec-relative-abundance} We begin by loading the relative abundance data and calculating dysbiosis scores for each sample based on Bray–Curtis dissimilarity from a healthy reference set. @@ -206,7 +207,7 @@ multimedia function to estimate path models for the treatment-to-mediator (α path), mediator-to-outcome (β path), and the combined indirect effect (α×β path). -We now fit the multimodal mediation model and inspect the overall indirect and +We now fit the multivariate mediation model and inspect the overall indirect and direct effects to understand whether microbiome composition mediates the treatment's effect on dysbiosis. @@ -225,7 +226,7 @@ exper <- mediation_data( mediators = medi_idx ) -# Fit the multimodal mediation model and estimate effects +# Fit the multivariate mediation model and estimate effects mdl <- multimedia(exper) res <- estimate(mdl, exper) @@ -369,7 +370,7 @@ ie_pw_fast <- tibble( mediator = top_meds, direct_setting = levels(tse$treatment)[1], contrast = paste(levels(tse$treatment)[1], "-", levels(tse$treatment)[2]), - indirect_effect = summary_df$estimate[match(top_meds, summary_df$mediator)] + indirect_effect = summary_df$effect[match(top_meds, summary_df$mediator)] ) # Visualize the selected mediators across samples @@ -383,7 +384,7 @@ same setup as before, using the time point as the treatment, dysbiosis score as the outcome, but pathways abundances as the mediation variables instead of the relative species abundance. -## Performing multimodal mediation analysis for pathways abundances {#sec-pathways-abundances} +## Performing multivariate mediation analysis for pathways abundances {#sec-pathways-abundances} We then repeat the same minimal preprocessing on the iHMP pathway abundance table (paired visits 1 and 27 for pathways). @@ -455,7 +456,7 @@ rownames(se) <- rownames(se) |> make.names(unique = TRUE) ``` -Next, we continue to fit the multimodal mediation model, and extract the overall +Next, we continue to fit the multivariate mediation model, and extract the overall and mediation-specific indirect effects. ```{r} @@ -473,7 +474,7 @@ exper <- mediation_data( mediators = medi_idx ) -# Fit the multimodal mediation model and estimate effects +# Fit the multivariate mediation model and estimate effects mdl <- multimedia(exper) res <- estimate(mdl, exper) @@ -505,7 +506,7 @@ both the overall and mediator-specific pathway indirect effects. set.seed(12345) # Bootstrap the overall indirect effect -boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 100) +boot_overall <- bootstrap(mdl, exper, c(indirect = indirect_overall), B = 10) # Summarize the bootstrap distribution and percentile interval summary(boot_overall$indirect) @@ -514,7 +515,7 @@ quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) # Plot the bootstrap distribution ggplot(boot_overall$indirect, aes(indirect_effect)) + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + - ggtitle("Bootstrap Distribution of Overall Indirect Effect (B = 100)") + + ggtitle("Bootstrap Distribution of Overall Indirect Effect (B = 10)") + labs(x = "Overall Indirect Effect", y = "Frequency") + theme_classic() ``` @@ -529,7 +530,7 @@ pathway abundance as well, with the point estimates and 95% CIs. set.seed(12345) # Bootstrap mediator-specific indirect effects -boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) +boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 10) # Remove bootstrap index boot_each <- boot_each$indirect[, -1] From 18a9c693850278808749a1ff178fceebf738f9c7 Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Sun, 9 Aug 2026 13:11:51 +0800 Subject: [PATCH 15/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 6c2b98e8..39fb086e 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -253,9 +253,13 @@ top_mediators <- effects_by_mediator |> print(top_mediators) ``` -To quantify uncertainty around the indirect effects, we use non-parametric -bootstrap resampling to compute **95% confidence intervals** for both the -overall and mediator-specific indirect effects. +To quantify uncertainty around the indirect effects, we use non-parametric +bootstrap resampling to compute **95% confidence intervals** for both the +overall and mediator-specific indirect effects. To keep the book build time +short, we use only `B = 10` bootstrap iterations below; for real applications +we recommend `B >= 500` (ideally 1000 or more) to obtain stable percentile +confidence intervals, so the CI widths shown here should be interpreted as +illustrative rather than definitive. ```{r} #| label: visualization_for_species_overall_indirect_effects @@ -273,12 +277,12 @@ quantile(boot_overall$indirect$indirect_effect, probs = c(0.025, 0.975)) # Plot the bootstrap distribution ggplot(boot_overall$indirect, aes(indirect_effect)) + geom_histogram(bins = 20, fill = "#69b3a2", color = "black") + - ggtitle("Bootstrap Distribution of Overall Indirect Effect (B = 100)") + + ggtitle("Bootstrap Distribution of Overall Indirect Effect (B = 10)") + labs(x = "Overall Indirect Effect", y = "Frequency") + theme_classic() ``` -Finally, we visualize the **mediator-specific indirect effects** using a forest +Finally, we visualize the **mediator-specific indirect effects** using a forest plot, displaying point estimates and 95% CIs to highlight the most influential mediators. @@ -301,7 +305,7 @@ indirect_each <- function(mdl, exper) { set.seed(12345) # Bootstrap mediator-specific indirect effects -boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 100) +boot_each <- bootstrap(mdl, exper, c(indirect = indirect_each), B = 10) # Remove bootstrap index boot_each <- boot_each$indirect[, -1] From 8aaf9f8ce336bef37856a2e45d3a252620f956ba Mon Sep 17 00:00:00 2001 From: YihanLiu4023 <149611140+YihanLiu4023@users.noreply.github.com> Date: Sun, 9 Aug 2026 13:29:54 +0800 Subject: [PATCH 16/16] Update multimedia.qmd --- inst/pages/multimedia.qmd | 34 ++++++++++++++-------------------- 1 file changed, 14 insertions(+), 20 deletions(-) diff --git a/inst/pages/multimedia.qmd b/inst/pages/multimedia.qmd index 39fb086e..89994887 100644 --- a/inst/pages/multimedia.qmd +++ b/inst/pages/multimedia.qmd @@ -26,7 +26,7 @@ insights into **how exposures impact outcomes via microbial features and** ```{r} #| label: fig_multivariate_mediation -#| fig-cap: Directed acyclic graph illustrating multivariate mediation, where an exposure affects many parallel mediators (here, illustrated as microbial features) within a single omic layer, which in turn affect the outcome. A direct path from exposure to outcome is also included (dashed). For tractability in high-dimensional mediation analysis, mediators are assumed to be conditionally independent given the exposure; this is a simplification, since real microbial features are often correlated. +#| fig-cap: Directed acyclic graph illustrating multivariate mediation where an exposure affects multiple mediators across modalities (microbial species and metabolites), which in turn affect the outcome. A direct path from exposure to outcome is also included (dashed). While gut microbial species biologically influence metabolite levels, for tractability in high-dimensional mediation analysis we assume species and metabolites act as conditionally independent parallel mediators given exposure, acknowledging this as an approximation #| echo: false library(ggdag) @@ -35,24 +35,24 @@ library(ggplot2) # Define DAG: exposure -> 4 parallel mediators -> outcome, plus a direct path dag <- dagify( - outcome ~ m1 + m2 + m3 + m4 + exposure, - m1 ~ exposure, - m2 ~ exposure, - m3 ~ exposure, - m4 ~ exposure, + outcome ~ sp1 + sp2 + met1 + met2 + exposure, + sp1 ~ exposure, + sp2 ~ exposure, + met1 ~ exposure, + met2 ~ exposure, exposure = "exposure", outcome = "outcome", labels = c( exposure = "Exposure", outcome = "Outcome", - m1 = "Feature 1", - m2 = "Feature 2", - m3 = "Feature 3", - m4 = "Feature 4" + sp1 = "Species 1", + sp2 = "Species 2", + met1 = "Metabolite 1", + met2 = "Metabolite 2" ), coords = list( - x = c(exposure = 0, m1 = 1, m2 = 1, m3 = 1, m4 = 1, outcome = 2), - y = c(exposure = 0, m1 = 1.5, m2 = 0.5, m3 = -0.5, m4 = -1.5, outcome = 0) + x = c(exposure = 0, sp1 = 1, sp2 = 1, met1 = 1, met2 = 1, outcome = 2), + y = c(exposure = 0, sp1 = 1.5, sp2 = 0.5, met1 = -0.5, met2 = -1.5, outcome = 0) ) ) @@ -215,15 +215,12 @@ treatment's effect on dysbiosis. #| label: overall_mediation_analysis #| message: false -# Use indices for mediators -medi_idx <- seq_len(nrow(tse)) - # Create the mediation data object exper <- mediation_data( tse, outcomes = "dysbiosis", treatments = "treatment", - mediators = medi_idx + mediators = rownames(tse) ) # Fit the multivariate mediation model and estimate effects @@ -467,15 +464,12 @@ and mediation-specific indirect effects. #| label: mediation_analysis #| message: false -# Use indices for mediators -medi_idx <- seq_len(nrow(se)) - # Create the mediation data object exper <- mediation_data( se, outcomes = "dysbiosis", treatments = "treatment", - mediators = medi_idx + mediators = rownames(se) ) # Fit the multivariate mediation model and estimate effects