-
Notifications
You must be signed in to change notification settings - Fork 44
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
maaslin2 --> maaslin3 #656
base: devel
Are you sure you want to change the base?
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -1136,14 +1136,17 @@ @Article{Das2020 | |
microbiomes with bias .pdf} | ||
} | ||
|
||
@Manual{Mallick2020, | ||
title = {MaAsLin 2: Multivariable Association in | ||
Population-scale Meta-omics Studies.}, | ||
author = {Himel Mallick and Ali Rahnavard and Lauren | ||
J. McIver}, | ||
year = {2020}, | ||
@Manual{Mallick2024, | ||
title = {MaAsLin 2: Refining and extending generalized | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This says maaslin 2 |
||
multivariable linear models for meta-omic association | ||
discovery.}, | ||
author = {Nickols, William A. and Kuntz, Thomas and Shen, | ||
Jiaxian and Maharjan, Sagun and Mallick, | ||
Himel and Franzosa, Eric A. and Thompson, | ||
Kelsey N. and Nearing, Jacob T. and Huttenhower, Curtis}, | ||
year = {2024}, | ||
note = {R/Bioconductor package}, | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It is not Bioconductor package |
||
url = {http://huttenhower.sph.harvard.edu/maaslin2}, | ||
url = {https://huttenhower.sph.harvard.edu/maaslin3/}, | ||
} | ||
|
||
@Article{Paulson2017, | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -316,50 +316,67 @@ ancombc2_out$res |> | |
kable() | ||
``` | ||
|
||
### MaAsLin2 | ||
### MaAsLin3 | ||
|
||
Let us next illustrate MaAsLin2 [@Mallick2020]. This method is based on | ||
Let us next illustrate MaAsLin3 [@Mallick2024]. This method is based on | ||
generalized linear models and flexible for different study designs | ||
and covariate structures. For details, check their | ||
[Biobakery tutorial](https://github.com/biobakery/biobakery/wiki/maaslin2). | ||
[Biobakery tutorial](https://github.com/biobakery/biobakery/wiki/MaAsLin3). | ||
|
||
```{r run-maaslin2, warning=FALSE, results="hide"} | ||
```{r run-maaslin3, warning=FALSE, results="hide"} | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. When you modify these chunk options, please make them to use quarto syntax. --> that way we gradually adopt the full quarto syntax in the book |
||
# Load package | ||
library(Maaslin2) | ||
library(maaslin3) | ||
|
||
# maaslin expects features as columns and samples as rows | ||
# for both the abundance table as well as metadata | ||
# MaAsLin3 takes tse as input data | ||
|
||
# We can specify different GLMs/normalizations/transforms. | ||
# specifying a ref is especially important if you have more than 2 levels | ||
maaslin2_out <- Maaslin2( | ||
input_data = as.data.frame(t(assay(tse))), | ||
input_metadata = as.data.frame(colData(tse)), | ||
output = "DAA example", | ||
transform = "AST", | ||
fixed_effects = "patient_status", | ||
# you can also fit MLM by specifying random effects | ||
# random_effects = c(...), | ||
reference = "patient_status,Control", | ||
normalization = "TSS", | ||
standardize = FALSE, | ||
# filtering was previously performed | ||
min_prevalence = 0) | ||
maaslin3_out <- maaslin3(input_data = tse, | ||
output = "DAA_example", | ||
formula = '~ patient_status', | ||
normalization = 'TSS', | ||
transform = 'LOG', | ||
augment = TRUE, | ||
standardize = TRUE, | ||
max_significance = 0.1, | ||
median_comparison_abundance = TRUE, | ||
median_comparison_prevalence = FALSE, | ||
max_pngs = 10, | ||
cores = 1, | ||
Comment on lines
+335
to
+344
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Use either " or ', not both. I think " is more common in the book so use it
Comment on lines
-343
to
+344
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. maaslin3 seems to have lots of different options. However, here we just want minimum example. That would be probably:
|
||
save_models = TRUE, | ||
verbosity = "WARN") | ||
``` | ||
|
||
Which genera are identified as differentially abundant? | ||
(leave out "head" to see all). | ||
|
||
```{r, maaslin2-res} | ||
maaslin2_out$results |> | ||
filter(qval < 0.05) |> | ||
kable() | ||
```{r, maaslin3-res} | ||
signif_results <- read.csv('DAA_example/significant_results.tsv', | ||
sep='\t') | ||
Comment on lines
+353
to
+354
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This does not look good. The indentation is not good There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Is this file the only format for results? I tried to quickly check, but it seems that the summary table is only available in the file |
||
head(signif_results, 20) %>% | ||
dplyr::mutate_if(is.numeric, .funs = function(x){( | ||
as.character(signif(x, 3)))}) %>% | ||
knitr::kable() %>% | ||
kableExtra::kable_styling("striped", position = 'center') %>% | ||
kableExtra::scroll_box(width = "800px", height = "400px") | ||
Comment on lines
+355
to
+360
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think this is too complex and ordinary reader might not understand what is happening here. I think you could just print the table with kable There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Moreover, elsewhere we use functions without defining the package, i.e., knitr::kable() --> kable() |
||
``` | ||
|
||
This will create a folder that is called like in the output specified | ||
This will create a folder that is called in the output specified | ||
above. It contains also figures to visualize difference between | ||
significant taxa. | ||
|
||
```{r, summary-plot} | ||
knitr::include_graphics("DAA_example/figures/summary_plot.png") | ||
Comment on lines
+367
to
+368
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This is ok to plot |
||
``` | ||
|
||
```{r, echo=FALSE, fig.show='hold', cache = FALSE, include=FALSE, eval=FALSE} | ||
prefix <- "DAA_example/figures/association_plots/patient_status" | ||
plot_vec <- c("/linear/patient_status_Coprobacter_linear.png", | ||
"/linear/patient_status_Faecalibacterium_linear.png", | ||
"/logistic/patient_status_Erysipelatoclostridium_logistic.png", | ||
) | ||
Comment on lines
+373
to
+376
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This does not run There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Error in c("/linear/patient_status_Coprobacter_linear.png", "/linear/patient_status_Faecalibacterium_linear.png", : There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Moreover, I think we could skip plotting these additional Maaslin3 plots. |
||
knitr::include_graphics(paste0(prefix, plot_vec)) | ||
``` | ||
|
||
### PhILR | ||
|
||
PhILR is a tree-based method that tests group-wise associations based on | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -34,7 +34,7 @@ gsEasy | |
igraph | ||
IntegratedLearner | ||
knitr | ||
Maaslin2 | ||
maaslin3 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This should be also added to DESCRIPTION to remotes |
||
mia | ||
miaTime | ||
miaViz | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Moreover, cite this paper instead
https://www.biorxiv.org/content/10.1101/2024.12.13.628459v1