Create multimedia.qmd - #768
Conversation
There was a problem hiding this comment.
Thanks!
Check the suggestions and the following:
-
increase support for the mia and (Tree)SummarizedExperiment methods when possible, this is very useful for compatibility with many other methods; I added suggestions. Also multimedia seems to support SE so this should be done.
-
Let us try to avoid repetition in the code; this has now two almost identical parts. Shall we prepare example data that is readily usable for the examples?
-
I am wondering if there would be a way to combine the two mediation chapters into one?
-
If we consider these methods widely usable we should consider providing wrappers for the interpretation, summaries, and visualizations of the results; either as contributions to "multimedia" package, or into mia or some other pkg (doesn't matter which pkg as long as it works)
|
Could you resolve the cases you have closed (from the "Resolve" buttons)? |
|
Fixed most of the detailed comments, still working on the more genral ones like creating pkg. Will finish this asap |
|
Can you resolve the completed suggestions above and confirm if this is ready to merge? |
Just fixed some minor inconsistence in the name of variables
|
To be checked in more detail before merge, regarding the following:
|
Hi Leo, thank you for the suggestions. I have added some wrapper functions, especially to where you indicated. For the data cleanning steps, do you mean we can remove those part, and add some texts to explain instead? |
|
The data cleaning part would work best if we include a suitable cleaned demo data in the R pkg, then we can skip the data cleaning steps in OMA. The problem is that if we could do this with every single method but then more than half of the book might easily be data cleaning examples, and this is also shifting the focus from the actual method to general data processing steps -> can we include demo data set/s in the multimedia, mia, or other pkg, or can we use already existing demo data sets? |
|
Regarding the wrappers, we could see if we can include these in a package if they would be generally useful anyway. |
|
Hi Leo – just wanted to quickly chime in here. This dataset is a special case where we don’t have a conventional, well-defined outcome for mediation analysis. Instead, the outcome is a dysbiosis score derived from taxonomic profiles, as described by Lloyd-Price et al. (2019), which, to my knowledge, hasn’t been implemented elsewhere. We are already drawing the data from curatedMetagenomicData, but creating the dysbiosis score is a necessary step to showcase the mediation analysis. |
|
Thanks. I think that dysbiosis score could then stay there. Let's check @TuomasBorman feedback. |
TuomasBorman
left a comment
There was a problem hiding this comment.
Thanks, looks very good! However, there are couple points to discuss.
This book is intended as teaching material to demonstrate specific analyses and tasks. The main focus is on the underlying ideas and concepts rather than implementation details. Therefore, the code is kept as simple as possible to achieve the intended outcomes. More polished analyses, including advanced plots, are better suited for workflow packages (we are also considering creating workflow package to showcase more complex analyses).
There are couple lengthy code chunks that do operations to results. Having these long chunks have opposite effect as readers are exhausted by the amount of code and they cannot focus on the main points. Everything that is not relevant for the main point should be removed even though plots or names would not be optimal.
Ideally, these wrappers will be implemented in a package, but of course it takes some time. At least, we should have plan to implement them. If we just leave it here, these will be never implemented as we are very busy with everything.
|
Can you @YihanLiu4023 confirm when you're ready with the updates? Can you also press the "Update branch" button above to ensure that this PR is in sync with the latest devel branch. |
|
Updated file submitted! |
antagomir
left a comment
There was a problem hiding this comment.
Thanks - some points to clarify still.
Overall this looks good but it will be essential to consider the following:
-
Data cleaning operations should be either ignored, or readily cleaned demo data set should be used (it can be placed to one of our packages). Otherwise OMA would be primarily a data cleaning example collection..
-
We need to check if any ready-made functionality is available to fetch and visualize results from mediation analyses. Now this is done with custom code. Maybe multimedia package has some utilities and if not, we should see if some tasks are so central that they should be standardized (either into multimedia package as contributed PR, or into mia).
| # Run mediation (non-delta) # | ||
| ############################# | ||
|
|
||
| # Convert the TreeSummarizedExperiment to be SummarizedExperiment |
There was a problem hiding this comment.
Why is this needed?
Note that TreeSE is an instance of SE
There was a problem hiding this comment.
I just realized here is an issue:

It looks like the package is kind of hard-coded when recognizing the types - it does not handle TreeSummarizedExperiment in mediation_data() and cannot recognizes it as SummarizedExperimentin this situation. Thus it seems like we still need to convert to SummarizedExperiment before running the mediation.
There was a problem hiding this comment.
I simplied a little bit on this as the following:
se_relative <- as(tse, "SummarizedExperiment")
if (is.null(rownames(se_relative))) rownames(se_relative) <- rownames(rowData(se_relative)) <- make.names(rownames(tse), unique = TRUE)
There was a problem hiding this comment.
The multimedia developer kindly fixed this for us: krisrs1128/multimedia#7
Now you can directly input TreeSE to multimedia_data
| # 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 |
There was a problem hiding this comment.
Can we prepare demo data in one of the R packages (e.g. mia) that is sufficiently clean?
OMA is dedicated to show how methods work, and data cleaning is something that is so common that it will overwhelm all examples in the book if we are not actively avoiding it.
There was a problem hiding this comment.
Yes, I can do that. Do I just need to create a new pull request to the mia package for the cleaned demo data?
There was a problem hiding this comment.
Yes, a PR. It should follow similar conventions than the other demo data sets.
There was a problem hiding this comment.
ibdhdb is already a mia dataset so it is not ideal to create a second version of this dataset just for the multimedia chapter.
Instead, let's simplify and minimise preprocessing just enough so that multimedia doesn't complain.
@YihanLiu4023, is it necessary to clean the pathway names? If it is necessary, could we simply replace all dots and underscores with space or another accepted character, in a single line of code?
| 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") | ||
|
|
There was a problem hiding this comment.
Can we consider having default visualization method for mediation analyses, for isntance in miaViz pkg (if not multimedia pkg itself)?
There was a problem hiding this comment.
We can add a package-native mediator visualization using plot_mediators(), which shows the observed outcome–mediator relationships (colored by treatment) for the top pathwise indirect effects. Here is the example plots for species:

(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.)
While a miaViz-based visualization would instead emphasize the taxonomic/pathway abundance structure (e.g., distributions/compositional patterns) rather than mediation-effect estimates, our original bootstrap forest plot complements plot_mediators() by summarizing effect uncertainty via CIs across mediators. Do you want plot_mediators() to be the primary/default visualization (with the forest plot as optional), or should we keep both as complementary views?
There was a problem hiding this comment.
We could include both for now imo,
There was a problem hiding this comment.
No problem, I would keep both for now.
|
This seems nearly complete. I picked the following points from above. Can you confirm this is OK @YihanLiu4023 ?
|
|
For the data cleanning, I'm preparing the data demo and will pull request to the mia package soon. For the visualization part, I just upload a version that keeps both of the things we discussed above. I would also uplaod a cleaned version (using the updated data demo in mia) very soon after the PR for data is merged. Thank you! |
|
I have create a PR to the mia package for the data demo, but it seems like the data file is not named correctly (may showing "Add files via upload"), please let me know if you cannot open it. |
@0xMuluh could you have a look? |
|
To complete soon. |
|
Hi! After meeting @himelmallick at the office, I had a brief look at this chapter and it looks very nice and almost ready! @TuomasBorman @antagomir, from your side what still needs to be done to finalise the chapter? From my side, I would have a couple of suggestions to simplify the code in addition to the others' comments above.
Hope this helps, and let me know if you need more feedback! :) |
|
I agree about these updates.
2-3. Agreed, co-pilots can probably handle much of this semi-automatically. How shall we proceed, is @YihanLiu4023 able to finalise these, or support needed? |
|
For the records, this is the reference list of miaViz: https://microbiome.github.io/miaViz/reference/index.html For here, at least the vignettes of plotMediation, plotForest and plotHistogram might be relevant to check. |
|
Thanks. It is also possible that (some of) these might not be directly applicable but the general idea is to streamline the analyses by providing easy-to-use default tools for standard use cases. The current case is also helping to identify where gaps/needs may exist. |
|
Hi, I uploaded a new version according to the comments but had 1 failing check: Style code. I also ran the style workflow locally in R. However, the automated styler check reformatted the R chunks back to 2-space indentation. Is there any way to solve this? |
|
Thanks! Looks pretty good! I still have a couple of comments. |
| #| label: load_iHMP_relative_demo_data | ||
| #| message: false | ||
|
|
||
| load("/Users/lyh_ciel/Library/CloudStorage/OneDrive-YaleUniversity/WCM/OMA/multimedia_data_demo.rda") |
There was a problem hiding this comment.
Instead of importing a ready-made object, preprocessing steps should be kept here as previously, similar to the other OMA chapters.
There was a problem hiding this comment.
P.S., now it should be possible to pass TreeSE to multimedia without having to convert it to SE: krisrs1128/multimedia#7 (comment)
There was a problem hiding this comment.
Sorry for the confusion, I forgot to change the path here while testing. Actually, I uploaded the preprocessed data to the mia pacakge with a new pull request as discussed before. I'm just wondering if we still want to load it through mia or we can just keep the data preprocessing here as previously?
There was a problem hiding this comment.
No problem, sorry I wasn't aware of it.
You are using the ibdhdb dataset from HMP2, correct? If so, it's already in mia: https://microbiome.github.io/mia/reference/ibdmdb.html.
You could import it with data("ibdmdb", package = "mia") and preprocess it for the multimedia analysis.
There was a problem hiding this comment.
Linking to the related PR: microbiome/mia#802 (comment)
There was a problem hiding this comment.
Thank you for this update! That's exactly the PR I created for this data. But this one is already preprocessed, so I may load it directly. Please let me know if you think it's better to show some preprocessings here in multimedia, I can make some changes.
| mediators = medi_idx | ||
| ) | ||
|
|
||
| mdl <- multimedia(exper) |
There was a problem hiding this comment.
While the three-line headers were not needed, the actual comments describing each line of code should be kept here and elsewhere.
There was a problem hiding this comment.
I added back some comments for understanding, but please let me know if you think more comments need to be added
There was a problem hiding this comment.
Thank you! If possible, please add comments to all lines of code, there are still some parts without comments.
There was a problem hiding this comment.
Also, remove full points (.) at the end of the comments
There was a problem hiding this comment.
I've added as many comments as possible, while also trying to keep things concise
RiboRings
left a comment
There was a problem hiding this comment.
Thank you so much, I'm confident we'll be ready in one or two more rounds of review. Please keep up the great work! 🥇
| 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) | ||
| } |
There was a problem hiding this comment.
As mentioned above, is this preprocessing necessary? If it is, can we try to reduce it to the minimum, 1/2 lines?
There was a problem hiding this comment.
It is only used to make pathway labels more readable in the plot. If you think it's okay, I can delete them. But currently, I just reduced this to two inline cleanup lines and removed the helper function.
There was a problem hiding this comment.
Thanks! It's this one, right?
pwy <- gsub("PWY0\\.|PWY\\.", "PWY-", rownames(summary_df))
pwy <- trimws(gsub("\\s+", " ", gsub("[_\\.]+", " ", pwy)))
rownames(summary_df) <- pwy
summary_df$mediator <- pwy
How about using stringr for all string processing in the chapter? That way would be more readable and it also supports pipe operations. For example, the code above becomes:
summary_df$mediator <- rownames(summary_df) |>
str_replace_all("PWY0\\.|PWY\\.", "PWY-") |>
str_remove_all("[_\\.]+") |>
str_remove_all("\\s+") |>
str_trim()
There was a problem hiding this comment.
To keep it simple, you can also check if the plot still looks good enough even when skipping some of those operations.
| diss <- as.matrix( | ||
| getDissimilarity( | ||
| tse, | ||
| method = "bray", | ||
| na.rm = TRUE, | ||
| assay.type = "relabundance" | ||
| ) | ||
| ) |
There was a problem hiding this comment.
Tidy operations are easier to read. For example:
# Calculate Bray-Curtis dissimilarity across samples
dist <- tse |>
getDissimilarity(...) |>
as.matrix()
Could you try to make this chunk (load_iHMP_relative_demo_data) more user-friendly through more tidyverse/dplyr verbs and |> pipes when possible?
There was a problem hiding this comment.
I just uploaded a version that load dataset after the preprocessing from mia, if you want to include some data cleanning here, I will revise this part later.
There was a problem hiding this comment.
In my understanding, the dataset you are using is ibdhmp from HMP2, correct? If so, we should use the original dataset, and apply all preprocessing steps in the multimedia chapter (while also keeping them minimal).
Please use:
# Import dataset
data("ibdmdb", package = "mia")
mae <- ibdmdb
| mediators = medi_idx | ||
| ) | ||
|
|
||
| mdl <- multimedia(exper) |
There was a problem hiding this comment.
Also, remove full points (.) at the end of the comments
|
Sorry I just realized the data loading is not yet updated (all the other comments till now is updated). I would push a new one right away. |
| # Use syntactically valid mediator names in model formulas | ||
| raw <- sub(".*s__", "", rownames(tse)) | ||
| clean1 <- str_remove_all(raw, "\\[|\\]") | ||
| clean2 <- str_replace_all(clean1, "[: \\.,]", "_") |
There was a problem hiding this comment.
clean1 and clean2 are used to make the mediator names formula-safe and more readable. I can delete this if you think this cleanup is unnecessary.
There was a problem hiding this comment.
It's ok to keep this, but as I commented above, let's make tidy and more readable using stringr verbs and the pipe:
# Use syntactically valid mediator names in model formulas
rownames(tse) <- rownames(tse) |>
str_remove(".*s__") |>
str_remove_all("\\[|\\]") |>
str_replace_all("[: \\.,]", "_") |>
make.names(unique = TRUE)
Change accordingly here and elsewhere.
| #| 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) |
There was a problem hiding this comment.
No need if it's too complicated, but would it be possible to plot this DAG with the ggdag package instead of DiagrammeR? That would align it with the OMA chapter on mediation.
|
|
||
| # Select the top mediators for the default mediator-level plot | ||
| top_meds <- summary_df %>% | ||
| dplyr::arrange(dplyr::desc(abs(estimate))) %>% |
There was a problem hiding this comment.
Remove package:: from function calls across the chapter unless necessary.
|
|
||
| 7. Repeat the process for the iHMP microbial pathways. | ||
|
|
||
| ### Performing multimodal mediation analysis for iHMP species relative abundance {#sec-relative-abundance} |
There was a problem hiding this comment.
Shouldn't this be a second-level header (##)? Same also for the next section.
| pwy <- gsub("PWY0\\.|PWY\\.", "PWY-", rownames(summary_df)) | ||
| pwy <- trimws(gsub("\\s+", " ", gsub("[_\\.]+", " ", pwy))) | ||
| rownames(summary_df) <- pwy | ||
| summary_df$mediator <- pwy |
There was a problem hiding this comment.
id.var argument of plotForest is "rownames" by default, so I don't think you need this line (457).
|
@YihanLiu4023 can you resolve the remaining comments, then we approach completion? |
|
@YihanLiu4023 kindly confirm when this is complete based on the feedback. |
No description provided.