-
Notifications
You must be signed in to change notification settings - Fork 72
MaAsLin3
MaAsLin3 is the next generation of MaAsLin (Microbiome Multivariable Associations with Linear Models). This comprehensive R package efficiently determines multivariable associations between clinical metadata and microbial meta-omics features. Relative to MaAsLin 2, MaAsLin 3 introduces the ability to quantify and test for both abundance and prevalence associations while accounting for compositionality. By incorporating generalized linear models, MaAsLin 3 accomodates most modern epidemiological study designs including cross-sectional and longitudinal studies.
If you use the MaAsLin 3 software, please cite our manuscript:
William A. Nickols, Jacob T. Nearing, Kelsey N. Thompson, Jiaxian Shen, Curtis Huttenhower. MaAsLin 3: Refining and extending generalized multivariate linear models for meta-omic association discovery. (In progress).
The user manual can be found here. If you have questions, please direct them to the MaAsLin 3 Forum.
- 1. Installing R
- 2. Installing MaAsLin 3
- 3. Microbiome association detection with MaAsLin 3
- 4. Advanced Topics
- 5. Command line
R is a programming language specializing in statistical computing and graphics. You can use R just the same as any other programming languages, but it is most useful for statistical analyses, with well-established packages for common tasks such as linear modeling, 'omic data analysis, machine learning, and visualization. We have created a short introduction tutorial here, but we note there are many other more comprehensive tutorials out there to learn R.
You can download and install the free R software environment here. Note that you should download the latest release - this will ensure the R version you have is compatible with MaAsLin 3.
RStudio is a freely available IDE (integrated development environment) for R. It is a "wrapper" around R with some additional functionalities that make programming in R a bit easier. Feel free to download RStudio and explore its interface - but it is not required for this tutorial.
If you already have R installed, then it is possible that the R version you have does not support MaAsLin 3. The easiest way to check this is to launch R/RStudio, and in the console ("Console" pane in RStudio), type in the following command (without the >
symbol):
sessionInfo()
The first line of output message should indicate your current R version. For example:
> sessionInfo()
R version 4.4.1 (2024-06-14)
For MaAsLin 3 to install, you will need R >= 4.3. If your version is older than that, please refer to section Installing R for the first time to download the latest R. Note that RStudio users will need to have RStudio "switch" to this later version once it is installed. This should happen automatically for Windows and Mac OS users when you relaunch RStudio. For Linux users you might need to bind the correct R executable. For more information, look here. Either way, once you have the correct version installed, launch the software and use sessionInfo()
to make sure that you indeed have R >= 4.3.
The latest development version of MaAsLin 3 can be installed from GitHub using the devtools
package.
# Install devtools if not present
if (!require('devtools', character.only = TRUE)) {
install.packages('devtools')
}
# Install MaAsLin 3
library("devtools")
install_github("biobakery/MaAsLin3")
for (lib in c('maaslin3', 'dplyr', 'ggplot2', 'knitr', 'kableExtra')) {
suppressPackageStartupMessages(require(lib, character.only = TRUE))
}
To run MaAsLin 3, the user must supply a table of per-sample feature abundances (with zeros still included), a table of per-sample metadata, and a model specifying how the metadata should relate to the feature prevalence (how likely the feature is to be present or absent) and abundance (how much of the feature is there if it's there). MaAsLin 3 will return a table of associations including an effect size and p-value for each feature-metadatum association and a folder of visuals including a summary plot and diagnostic plots for significant associations.
MaAsLin3 requires two input files.
- Feature abundance data frame
- Formatted with features as columns and samples as rows.
- The transpose of this format is also okay.
- Possible features include taxonomy or genes. These can be relative abundances or counts.
- This can be a filepath to a tab-delimited file.
- Metadata data frame
- Formatted with variables as columns and samples as rows.
- The transpose of this format is also okay.
- Possible metadata include gender or age.
- This can be a filepath to a tab-delimited file.
The data file can contain samples not included in the metadata file (along with the reverse case). For both cases, those samples not included in both files will be removed from the analysis. Also, the samples do not need to be in the same order in the two files.
Example input files can be found in the inst/extdata
folder of the MaAsLin 3 source. The files provided were generated from the Human Microbiome Project 2 (HMP2) data which can be downloaded from https://ibdmdb.org/.
-
HMP2_taxonomy.tsv
: a tab-delimited file with samples as rows and species as columns. It is a subset of the full HMP2 taxonomy that includes just some of the the species abundances. -
HMP2_metadata.tsv
: a tab-delimited file with samples as rows and metadata as columns. It is a subset of the full HMP2 metadata that includes just some of the fields.
# Read abundance table
taxa_table_name <- system.file("extdata", "HMP2_taxonomy.tsv", package = "maaslin3")
taxa_table <- read.csv(taxa_table_name, sep = '\t')
# Read metadata table
metadata_name <- system.file("extdata", "HMP2_metadata.tsv", package = "maaslin3")
metadata <- read.csv(metadata_name, sep = '\t')
# Factor the categorical variables to test IBD against healthy controls
metadata$diagnosis <-
factor(metadata$diagnosis, levels = c('nonIBD', 'UC', 'CD'))
metadata$dysbiosis_state <-
factor(metadata$dysbiosis_state, levels = c('none', 'dysbiosis_UC', 'dysbiosis_CD'))
metadata$antibiotics <-
factor(metadata$antibiotics, levels = c('No', 'Yes'))
# Check how the dataframes look like
taxa_table[1:5, 1:5]
metadata[1:5, 1:5]
Phocaeicola_vulgatus Faecalibacterium_prausnitzii Bacteroides_uniformis Prevotella_copri_clade_A Bacteroides_stercoris
CSM5FZ3N_P 0.4265226 0.060255109 0.2692411314 0.00000000000 0.000000000
CSM5FZ3R_P 0.5369584 0.007396904 0.2526048469 0.00000000000 0.008390958
CSM5FZ3T_P 0.5911821 0.000000000 0.0000000000 0.00000000000 0.000000000
CSM5FZ3V_P 0.2661378 0.029680329 0.4004265260 0.00000000000 0.000000000
CSM5FZ3X_P 0.6601039 0.003596740 0.0008804283 0.00001308102 0.001335669
participant_id site_name week_num reads diagnosis
CSM5FZ3N_P C3001 Cedars-Sinai 0 9961743 CD
CSM5FZ3R_P C3001 Cedars-Sinai 2 16456391 CD
CSM5FZ3T_P C3002 Cedars-Sinai 0 10511448 CD
CSM5FZ3V_P C3001 Cedars-Sinai 6 17808965 CD
CSM5FZ3X_P C3002 Cedars-Sinai 2 13160893 CD
MaAsLin 3 is run by specifying a param_list
that contains the abundance table (input_data
), the metadata table (input_metadata
), the output directory (output
), and a model. The model can come from a formula or vectors of terms. In either case, variable names should not have spaces or unusual characters.
-
Formula: The
formula
parameter should be set to any formula that satisfies thelme4
specifications: fixed effects, random effects (section 4.2), interaction terms (section 4.3), polynomial terms, and more can all be included. If categorical variables are included as fixed effects, each level will be tested against the first factor level. In addition, ordered predictors (section 4.4) and group predictors (section 4.5) can be included by includinggroup(variable_name)
andordered(variable_name)
in the formula. -
Vectors: Alternatively, a vector of variable names can be supplied to the parameters
fixed_effects
,random_effects
,group_effects
, andordered_effects
. Categorical variables should either be ordered as factors beforehand, orreference
should be provided as a string of 'variable,reference' semi-colon delimited for multiple variables (e.g.,variable_1,reference_1;variable_2,reference_2
). NOTE: adding a space between the variable and level might result in the wrong reference level being used.
Because MaAsLin 3 identifies prevalence (presence/absence) associations, untransformed read depth (number of reads in the sample) should be included as a covariate if available. Deeper sequencing will likely increase feature detection in a way that could spuriously correlate with metadata of interest when read depth is not included in the model. The read depth should be included untransformed (i.e., not log read depth) since the log odds of a microbe being missed by chance scale approximately linearly with the untransformed read depth.
The following command runs MaAsLin 3 on the HMP2 data, running a multivariable
regression model to test for the association between microbial species abundance and prevalence
versus IBD diagnosis and dysbiosis scores while controlling for antibiotic usage, age, and the sampling depth (reads). The abundance associations come from fitting multivariate linear models to the log base 2 relative abundances of features in samples in which the feature is present. The prevalence associations come from fitting multivariate logistic models to the presence/absence data for a feature across all samples. Outputs are directed to a folder called hmp2_output
under the current working directory (output = "hmp2_output"
).
To show one of each type of plot (see below), the maximum number of significant associations to plot has been increased with max_pngs=100
(default 30
). All other parameters beyond the first four are the default parameters but are still included for clarity:
- Total sum scaling (
normalization = 'TSS'
) with a log transformation ('transform = 'LOG'
) afterwards will be used. These are almost always the recommended choices (and all MaAsLin 3 evaluations were performed with these options), but other normalizations and transformations are allowed (seenormalization
andtransform
in?maaslin3
). - A data augmentation procedure is used to avoid linear separability in the logistic regressions (
augment = TRUE
). This is almost always recommended, though it can be turned off (seeaugment
in the manual). - Continuous metadatum variables are z-score standardized (subtract the mean, divide by the standard deviation) with
standardize = TRUE
so that the coefficients will be on the same scale (improving visualization). - A nominal FDR level of 0.1 (
max_significance = 0.1
) determines which associations are significant, but this can also be changed (seemax_significance
in the manual). - For the abundance coefficients, to account for compositionality in relative abundance data, significance is determined by comparing against the median coefficient for a metadatum across the features (
median_comparison_abundance
). Since prevalence coefficients do not have the same compositional properties, they are instead compared against 0 (median_comparison_prevalence
). See below for a discussion of these paremeters. - One CPU is used for fitting (
cores = 1
).
The abundance, prevalence, and variance filtering parameters are not included, and they are 0 by default. Different from other differential abundance tools, low prevalence features should not be filtered out since the prevalence modeling in MaAsLin 3 already accounts for high proportions of zeros.
param_list <- list(input_data = taxa_table,
input_metadata = metadata,
output = 'hmp2_output', # Maaslin3 creates the directory if it doesn't exist.
formula = '~ diagnosis + dysbiosis_state + antibiotics + age + reads',
normalization = 'TSS',
transform = 'LOG',
augment = TRUE,
standardize = TRUE,
max_significance = 0.1,
median_comparison_abundance = TRUE,
median_comparison_prevalence = FALSE,
max_pngs = 100,
cores = 1)
fit_out <- maaslin3(param_list)
When median_comparison_abundance
or median_comparison_prevalence
are on, the coefficients for a metadatum will be tested against the median coefficient for that metadatum (median across the features). Otherwise, the coefficients will be tested against 0. For abundance associations, this is designed to account for compositionality, the idea that testing against zero on the relative scale is not equivalent to testing against zero on the absolute scale. The user manual contains a more extensive discussion of when to use median_comparison_abundance
and median_comparison_prevalence
, but in summary:
-
median_comparison_abundance
isTRUE
by default and should be used to make inference on the absolute scale when using relative abundance data. Whenmedian_comparison_abundance
isTRUE
, only the p-values and q-values change; the coefficients returned are still the relative abundance coefficients. -
median_comparison_abundance
should beFALSE
when (1) testing for significant relative associations, (2) testing for absolute abundance associations under the assumption that the total absolute abundance is not changing, or (3) testing for significant absolute associations when supplying spike-in or total abundances withunscaled_abundance
. -
median_comparison_prevalence
isFALSE
by default.
The main output from MaAsLin 3 is the list of significant associations in significant_results.tsv
. This lists all associations that pass MaAsLin 3's significance
threshold, ordered by increasing joint q-value. The format is:
-
feature
andmetadata
are the feature and metadata names. -
value
andname
are the value of the metadata and variable name from the model. -
coef
andstderr
are the fit coefficient and standard error from the model. In abundance models, a one-unit change in the metadatum variable corresponds to a$2^{\textrm{coef}}$ fold change in the relative abundance of the feature. In prevalence models, a one-unit change in the metadatum variable corresponds to a$\textrm{coef}$ change in the log-odds of a feature being present. -
pval_individual
andqval_individual
are the p-value and false discovery rate (FDR) corrected q-value of the individual association. -
model
specifies whether the association is abundance or prevalence. -
N
andN.not.zero
are the total number of data points and the total number of non-zero data points for the feature. -
pval_joint
andqval_joint
are the p-value and q-value of the joint prevalence and abundance association. The p-value comes from plugging in the minimum of the association's abundance and prevalence p-values into the Beta(1,2) CDF. It is interpreted as the probability that either the abundance or prevalence association would be as extreme as observed if there was neither an abundance nor prevalence association between the feature and metadatum.
feature | metadata | value | name | coef | stderr | pval_individual | qval_individual | model | N | N.not.zero | pval_joint | qval_joint |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Escherichia_coli | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | 3.533 | 0.3868 | 0 | 0 | abundance | 1527 | 707 | 0 | 0 |
Veillonella_parvula | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | 4.055 | 0.4878 | 0 | 0 | abundance | 1527 | 547 | 0 | 0 |
Escherichia_coli | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | 1.583 | 0.2373 | 2.57e-11 | 2.921e-10 | prevalence | 1527 | 707 | 0 | 0 |
Veillonella_parvula | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | 0.8518 | 0.2146 | 0.00007222 | 0.0002447 | prevalence | 1527 | 547 | 0 | 0 |
Phocaeicola_sartorii | reads | reads | reads | -0.2403 | 0.1483 | 1 | 1 | abundance | 1527 | 424 | 1.083e-43 | 3.817e-41 |
Phocaeicola_sartorii | reads | reads | reads | 1.095 | 0.07872 | 5.417e-44 | 5.726e-41 | prevalence | 1527 | 424 | 1.083e-43 | 3.817e-41 |
Faecalibacterium_prausnitzii | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | -2.184 | 0.3095 | 2.913e-10 | 0.00000001638 | abundance | 1527 | 1374 | 7.086e-35 | 1.872e-32 |
Faecalibacterium_prausnitzii | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | -3.492 | 0.2822 | 3.543e-35 | 1.872e-32 | prevalence | 1527 | 1374 | 7.086e-35 | 1.872e-32 |
Bacteroides_uniformis | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | -3.123 | 0.3959 | 3.992e-13 | 3.992e-11 | abundance | 1527 | 1244 | 3.151e-31 | 6.662e-29 |
Bacteroides_uniformis | dysbiosis_state | dysbiosis_CD | dysbiosis_statedysbiosis_CD | -2.912 | 0.2493 | 1.576e-31 | 5.551e-29 | prevalence | 1527 | 1244 | 3.151e-31 | 6.662e-29 |
Bifidobacterium_longum | age | age | age | -0.1291 | 0.1063 | 1 | 1 | abundance | 1527 | 843 | 2.022e-29 | 3.562e-27 |
Bifidobacterium_longum | age | age | age | -0.6724 | 0.05939 | 1.011e-29 | 2.672e-27 | prevalence | 1527 | 843 | 2.022e-29 | 3.562e-27 |
MaAsLin 3 generates two types of output files explained below: data and visualizations. In addition, the object returned from maaslin3
(fit_out
above) contains all the data and results (see ?maaslin_fit
).
- Data output files
-
significant_results.tsv
- Described above.
-
all_results.tsv
- Same format as above but for all associations and with an extra column listing any errors from the model fitting.
-
features
- This folder includes the filtered, normalized, and transformed versions of the input feature table.
- These steps are performed sequentially in the above order.
- If an option is set such that a step does not change the data, the resulting table will still be output.
-
models_LM.rds
andmodels_logistic.rds
- These files contain a list with every model fit object (
LM
for linear models,logistic
for logistic models). - It will only be generated if
save_models
is set to TRUE.
- These files contain a list with every model fit object (
-
residuals_LM.rds
andresiduals_logstic.rds
- These files contain a data frame with residuals for each feature. For the logistic fits, the residuals are the deviance residuals.
-
fitted_LM.rds
andfitted_logistic.rds
- These files contain a data frame with fitted values for each feature.
-
ranef_LM.rds
andranef_logistic.rds
- These files contain a data frame with extracted random effects for each feature (when random effects are specified).
-
maaslin3.log
- This file contains all log information for the run.
- It includes all settings, warnings, errors, and steps run.
-
- Visualization output files
-
summary_plot.pdf
- This file contain a combined coefficient plot and heatmap of the most significant associations. In the heatmap, one star indicates the individual q-value is below the parameter
max_significance
, and two stars indicate the individual q-value is belowmax_significance
divided by 10.
- This file contain a combined coefficient plot and heatmap of the most significant associations. In the heatmap, one star indicates the individual q-value is below the parameter
-
-
association_plots/[metadatum]_[feature]_[association].pdf
- A plot is generated for each significant association up to
max_pngs
. - Scatter plots are used for continuous metadata abundance associations.
- Box plots are used for categorical data abundance associations.
- Box plots are used for continuous data prevalence associations.
- Grids are used for categorical data prevalence associations.
- Data points plotted are after filtering, normalization, and transformation so that the scale in the plot is the scale that was used in fitting.
- A plot is generated for each significant association up to
At the top right of each association plot is the name of the significant association in the results file, the FDR corrected q-value for the individual association, the number of samples in the dataset, and the number of samples with non-zero abundances for the feature. In the plots with categorical metadata variables, the reference category is on the left, and the significant q-values and coefficients in the top right are in the order of the values specified above. Because the displayed coefficients correspond to the full fit model with (possibly) scaled metadata variables, the marginal association plotted might not match the coefficient displayed. However, the plots are intended to provide an interpretable visual while usually agreeing with the full model.
There are a few common issues to check for in the results:
- When warnings or errors are thrown during the fitting process, they are recorded in the
error
column ofall_results.tsv
. Often, these indicate model fitting failures or poor fits that should not be trusted, but sometimes the warnings can be benign, and the model fit might still be reasonable. Users should check associations of interest if they produce errors. - Despite the error checking, significant results could still result from poor model fits. These can usually be diagnosed with the visuals in the
association_plots
directory.- Any significant abundance associations with a categorical variable should usually have at least 10 observations in each category.
- Significant prevalence associations with categorical variables should also have at least 10 samples in which the feature was present and at least 10 samples in which it was absent for each significant category.
- Significant abundance associations with continuous metadata should be checked visually for influential outliers.
- There are also a few rules of thumb to keep in mind.
- Models should ideally have about 10 times as many samples (all samples for logistic fits, non-zero samples for linear fits) as covariate terms (all continuous variables plus all categorical variable levels).
- Coefficients (effect sizes) larger than about 15 in absolute value are usually suspect unless very small unstandardized predictors are being included. (A coefficient of 15 corresponds to a fold change >30000!). If you encounter such coefficients, check that (1) no error was thrown, (2) the diagnostics look reasonable, (3) a sufficient number of samples were used in fitting, (4) the q-value is significant, (5) the metadata are not highly collinear, and (6) the random effects are plausible.
Once maaslin3
has been run once, maaslin_plot_results_from_output
can be run with the same param_list
as the original maaslin3
run to (re-)create the plots. This allows the user to plot the associations even without having the R object returned by maaslin_fit
or maaslin3
(e.g., if fitting models through the command line). It is recommended to fit models with simple variables names that are robust to formula plotting and then convert these into proper names for plotting. Likewise, heatmap_vars
and coef_plot_vars
can be specified in the param_list
at any point, but it is usually easier to see how the names come out first and then choose which metadata variables will go in the coefficient plot and which will go in the heatmap afterwards with maaslin_plot_results_from_output
.
# This section is necessary for updating the summary plot and the association plots
# Rename results file with clean titles
all_results <- read.csv('hmp2_output/all_results.tsv', sep='\t')
all_results <- all_results %>%
mutate(metadata = case_when(metadata == 'age' ~ 'Age',
metadata == 'antibiotics' ~ 'Abx',
metadata == 'diagnosis' ~ 'Diagnosis',
metadata == 'dysbiosis_state' ~ 'Dysbiosis',
metadata == 'reads' ~ 'Read depth'),
value = case_when(value == 'dysbiosis_CD' ~ 'CD',
value == 'dysbiosis_UC' ~ 'UC',
value == 'Yes' ~ 'Used', # Antibiotics
value == 'age' ~ 'Age',
value == 'reads' ~ 'Read depth',
TRUE ~ value),
feature = gsub('_', ' ', feature) %>%
gsub(pattern = 'sp ', replacement = 'sp. '))
# Write results
write.table(all_results, 'hmp2_output/all_results.tsv', sep='\t')
# Set the new heatmap and coefficient plot variables and order them
param_list$heatmap_vars = c('Dysbiosis UC', 'Diagnosis UC', 'Abx Used', 'Age', 'Read depth')
param_list$coef_plot_vars = c('Dysbiosis CD', 'Diagnosis CD')
# This section is necessary for updating the association plots
colnames(param_list$input_data) <- gsub('_', ' ', colnames(param_list$input_data)) %>%
gsub(pattern = 'sp ', replacement = 'sp. ')
# Rename the features in the norm transformed data file
filtered_data_norm_transformed <-
read.csv('hmp2_output/features/filtered_data_norm_transformed.tsv', sep='\t')
colnames(filtered_data_norm_transformed) <-
gsub('_', ' ', colnames(filtered_data_norm_transformed)) %>%
gsub(pattern = 'sp ', replacement = 'sp. ')
write.table(filtered_data_norm_transformed,
'hmp2_output/features/filtered_data_norm_transformed.tsv',
sep='\t', row.names = F)
# Rename the metadata like in the outputs table
colnames(param_list$input_metadata) <-
case_when(colnames(param_list$input_metadata) == 'age' ~ 'Age',
colnames(param_list$input_metadata) == 'antibiotics' ~ 'Abx',
colnames(param_list$input_metadata) == 'diagnosis' ~ 'Diagnosis',
colnames(param_list$input_metadata) == 'dysbiosis_state' ~ 'Dysbiosis',
colnames(param_list$input_metadata) == 'reads' ~ 'Read depth',
TRUE ~ colnames(param_list$input_metadata))
param_list$input_metadata <- param_list$input_metadata %>%
mutate(Dysbiosis = case_when(Dysbiosis == 'dysbiosis_UC' ~ 'UC',
Dysbiosis == 'dysbiosis_CD' ~ 'CD',
Dysbiosis == 'none' ~ 'None') %>%
factor(levels = c('None', 'UC', 'CD')),
Abx = case_when(Abx == 'Yes' ~ 'Used',
Abx == 'No' ~ 'Not used') %>%
factor(levels = c('Not used', 'Used')),
Diagnosis = case_when(Diagnosis == 'nonIBD' ~ 'non-IBD',
TRUE ~ Diagnosis) %>%
factor(levels = c('non-IBD', 'UC', 'CD')))
# Recreate the plots
scatter_plots <- maaslin_plot_results_from_output(param_list)
In the new summary plot below, we can see that the feature names are cleaned up, the metadata names are cleaned up, the set of metadata variables used in the coefficient plot is different, and the metadata used in the heatmap is reordered.
In the association plots, the taxa and metadata have been renamed to be consistent with the results file from earlier.
Most microbiome methodologies have historically focused on relative abundances (proportions out of 1). However, some experimental protocols can enable estimation of absolute abundances (cell count/concentration). MaAsLin 3 can be used with two types of absolute abundance estimation: spike-ins and total abundance scaling. In a spike-in procedure, a small, known quantity of a microbe that otherwise would not be present in the sample is added, and the sequencing procedure is carried out as usual. Then, the absolute abundance of a microbe already in the community is estimated as:
The following section shows a synthetic spike-in procedure with simulated data from SparseDOSSA2. The abundance table is like any other abundance table input to MaAsLin 3 except that 'Feature101' is the spike-in (which must be present in every sample). The scaling factor data frame has 'Feature101' as its only column name with the absolute abundance of the spike-in for each sample (rownames). If the same quantity of the spike-in was added to all samples, this column could be entirely 1 (or any other arbitrary value). When using a spike-in procedure, the scaling factor data frame should have a single column named the same as one of the features in the abundance table, which will be used as the spike-in.
# Abundance table
taxa_table_name <- system.file("extdata", "abundance_spike_in_ex.tsv", package = "maaslin3")
spike_in_taxa_table <- read.csv(taxa_table_name, sep = '\t')
# Metadata table
metadata_name <- system.file("extdata", "metadata_spike_in_ex.tsv", package = "maaslin3")
spike_in_metadata <- read.csv(metadata_name, sep = '\t')
for (col in c('Metadata_1', 'Metadata_2', 'Metadata_5')) {
spike_in_metadata[,col] <- factor(spike_in_metadata[,col])
}
# Spike-in table
unscaled_name <- system.file("extdata", "scaling_factors_spike_in_ex.tsv", package = "maaslin3")
spike_in_unscaled <- read.csv(unscaled_name, sep = '\t')
spike_in_taxa_table[c(1:5, 101),1:5]
spike_in_metadata[1:5,]
spike_in_unscaled[1:5, , drop=F]
Sample1 Sample2 Sample3 Sample4 Sample5
Feature1 0 0 0 0 0
Feature2 0 0 0 0 0
Feature3 554 0 1737 0 0
Feature4 0 0 0 0 0
Feature5 0 0 0 0 0
Feature101 15464 4788 5155 2052 11540
Metadata_1 Metadata_2 Metadata_3 Metadata_4 Metadata_5 ID
Sample1 1 1 -0.41011453 0.6352468 0 Subject1
Sample2 0 1 -1.35105918 0.2362369 0 Subject2
Sample3 1 1 -0.48798215 -0.5458514 0 Subject3
Sample4 0 0 -0.07268457 -0.9903428 1 Subject4
Sample5 1 0 0.27201198 1.1333223 1 Subject5
Feature101
Sample1 418
Sample2 183
Sample3 38
Sample4 3230
Sample5 243
The following code fits the absolute abundance model. Note that the normalization must be set to TSS
since the procedure involves scaling the relative abundance ratios. Also, median_comparison_abundance
is now set to FALSE
since we want to test the absolute coefficients against zero. The data frame of absolute abundances is included as the unscaled_abundance
parameter, and the spike-in strategy will automatically be detected based on the column name.
param_list <- list(input_data = spike_in_taxa_table,
input_metadata = spike_in_metadata,
output = 'spike_in_demo',
formula = '~ Metadata_1 + Metadata_2 + Metadata_3 + Metadata_4 + Metadata_5',
normalization = 'TSS',
transform = 'LOG',
median_comparison_abundance = FALSE,
unscaled_abundance = spike_in_unscaled)
fit_out <- maaslin3(param_list)
The absolute abundance coefficients can be inspected:
feature | metadata | value | name | coef | stderr | pval_individual | qval_individual | model | N | N.not.zero | pval_joint | qval_joint |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Feature58 | Metadata_2 | 1 | Metadata_21 | -4.897 | 0.5194 | 7.625e-13 | 2.28e-10 | abundance | 100 | 58 | 1.525e-12 | 6.786e-10 |
Feature58 | Metadata_2 | 1 | Metadata_21 | -0.6161 | 0.4647 | 0.1849 | 0.9788 | prevalence | 100 | 58 | 1.525e-12 | 6.786e-10 |
Feature88 | Metadata_5 | 1 | Metadata_51 | 7.86 | 0.8193 | 2.994e-12 | 4.475e-10 | abundance | 100 | 49 | 5.987e-12 | 0.000000001332 |
Feature88 | Metadata_5 | 1 | Metadata_51 | -0.3644 | 0.6519 | 0.5761 | 0.9991 | prevalence | 100 | 49 | 5.987e-12 | 0.000000001332 |
Feature3 | Metadata_1 | 1 | Metadata_11 | -1.479 | 0.9172 | 0.1162 | 0.7639 | abundance | 100 | 40 | 0.0000001174 | 0.00001742 |
Feature3 | Metadata_1 | 1 | Metadata_11 | 4.989 | 0.9201 | 0.00000005871 | 0.00002613 | prevalence | 100 | 40 | 0.0000001174 | 0.00001742 |
The following section shows a synthetic total abundance scaling procedure with simulated data from SparseDOSSA2. The abundance table is like any other abundance table input to MaAsLin 3 without any extra features. The scaling factor data frame has 'total' as its only column name with the total absolute abundance for each sample (rownames).
# Abundance table
taxa_table_name <- system.file("extdata", "abundance_total_ex.tsv", package = "maaslin3")
total_scaling_taxa_table <- read.csv(taxa_table_name, sep = '\t')
# Metadata table
metadata_name <- system.file("extdata", "metadata_total_ex.tsv", package = "maaslin3")
total_scaling_metadata <- read.csv(metadata_name, sep = '\t')
for (col in c('Metadata_1', 'Metadata_3', 'Metadata_5')) {
spike_in_metadata[,col] <- factor(spike_in_metadata[,col])
}
# Total abundance table
unscaled_name <- system.file("extdata", "scaling_factors_total_ex.tsv", package = "maaslin3")
total_scaling_unscaled <- read.csv(unscaled_name, sep = '\t')
total_scaling_taxa_table[1:5, 1:5]
total_scaling_metadata[1:5,]
total_scaling_unscaled[1:5, , drop=F]
Sample1 Sample2 Sample3 Sample4 Sample5
Feature1 11837 0 0 8 68
Feature2 0 0 0 0 1
Feature3 0 0 0 0 0
Feature4 0 0 1 0 0
Feature5 0 0 0 0 0
Metadata_1 Metadata_2 Metadata_3 Metadata_4 Metadata_5 ID
Sample1 0 1.2410844 0 0.4341439 1 Subject1
Sample2 1 0.0176896 1 1.4894610 1 Subject2
Sample3 1 -1.7370341 1 -0.7973592 1 Subject3
Sample4 1 0.1500793 0 3.2816068 0 Subject4
Sample5 0 -0.2819085 0 -0.6731906 0 Subject5
total
Sample1 1352.276
Sample2 1720.637
Sample3 64194.918
Sample4 50256.129
Sample5 1132.367
The following code fits the absolute abundance model. As before, TSS
must be set to TRUE
, median_comparison_abundance
is set to FALSE
, and the data frame of absolute abundances is included as the unscaled_abundance
parameter. The total abundance scaling procedure will be detected based on the column name.
param_list <- list(input_data = total_scaling_taxa_table,
input_metadata = total_scaling_metadata,
output = 'total_scaling_demo',
formula = '~ Metadata_1 + Metadata_2 + Metadata_3 + Metadata_4 + Metadata_5',
normalization = 'TSS',
transform = 'LOG',
median_comparison_abundance = FALSE,
unscaled_abundance = total_scaling_unscaled)
fit_out <- maaslin3(param_list)
The absolute abundance coefficients can be inspected:
feature | metadata | value | name | coef | stderr | pval_individual | qval_individual | model | N | N.not.zero | pval_joint | qval_joint |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Feature94 | Metadata_1 | Metadata_1 | Metadata_1 | 2.831 | 0.3736 | 3.322e-11 | 0.000000009568 | abundance | 100 | 94 | 6.644e-11 | 0.00000003056 |
Feature94 | Metadata_1 | Metadata_1 | Metadata_1 | -0.5294 | 0.4113 | 0.198 | 0.9853 | prevalence | 100 | 94 | 6.644e-11 | 0.00000003056 |
Feature80 | Metadata_3 | Metadata_3 | Metadata_3 | 2.86 | 0.3161 | 3.318e-10 | 0.00000004778 | abundance | 100 | 37 | 6.636e-10 | 0.0000001526 |
Feature80 | Metadata_3 | Metadata_3 | Metadata_3 | 0.4915 | 0.224 | 0.02822 | 0.5522 | prevalence | 100 | 37 | 6.636e-10 | 0.0000001526 |
Feature80 | Metadata_5 | Metadata_5 | Metadata_5 | 2.218 | 0.2961 | 0.00000001924 | 0.000001847 | abundance | 100 | 37 | 0.00000003847 | 0.000005899 |
Feature80 | Metadata_5 | Metadata_5 | Metadata_5 | 0.1294 | 0.2179 | 0.5526 | 0.9984 | prevalence | 100 | 37 | 0.00000003847 | 0.000005899 |
Certain studies have a natural "grouping" of sample observations, such as by subject in longitudinal designs or by family in family designs. It is important for statistic analysis to address the non-independence between samples belonging to the same group. MaAsLin 3 provides a simple interface for this with random effects, where the user can specify the grouping variable to run a mixed effect model. This grouping variable can be provided with the random_effects
parameter or specified in the model with (1|grouping_variable)
. This will add in a "random intercept," a per-group adjustment to the model intercept. More complicated random effects can be specified in the formula in accordance with lme4
formula parsing. For example, HMP2 has a longitudinal design where the same subject (column participant_id
) can have multiple samples. We thus use participant_id
as a random effect grouping variable in the model.
param_list <- list(input_data = taxa_table,
input_metadata = metadata,
output = 'random_effects_output',
formula = '~ diagnosis + dysbiosis_state + antibiotics + age + reads + (1|participant_id)',
plot_summary_plot = F,
plot_associations = F)
fit_out <- maaslin3(param_list)
Note that no participant_id
terms are included in the outputs; the random intercepts are just used to control for grouping. If you are interested in testing the effect of time in a longitudinal study, the time point variable should be included as a fixed effect in your MaAsLin 3 formula.
Mixed effects models take longer to fit because of their complexity, and the fitting may fail more often because of the extra terms. However, it is important to account for this non-independence, and MaAsLin 3 is still able to fit complex models with thousands of samples and thousands of features in a few hours (and increasing the CPUs used with cores
can help). However, random effects are not the solution to every correlated data problem. In particular, random intercept models may produce poor model fits when there are fewer than 5 observations per group. In these scenarios, per-group fixed effects should be used. For example, in a before-and-after treatment study with one sample before and one sample after, the participant IDs and a pre/post indicator should be included as fixed effects, but only the pre/post indicator should be retained for analysis afterwards.
One of the benefits of MaAsLin 3's formula mode is the ability to include interaction terms. Mathematically, an interaction term corresponds to the product of two columns in the design matrix. When a continuous variable is interacted with a categorical term, the interaction term corresponds to the change in the continuous variable's slope between the categories. For two categorical variables interacted, see below; it is better explained through an example. In the formula, interactions can be specified with the :
symbol to include only the interaction term(s) or the *
symbol to include both the interaction term and the non-interacted terms (i.e., a*b
adds the terms a
, b
, and a:b
).
param_list <- list(input_data = taxa_table,
input_metadata = metadata,
output = 'interaction_output',
formula = '~ sex*diagnosis + antibiotics + age + reads')
fit_out <- maaslin3(param_list)
In the model above, we have specified an interaction between sex
and diagnosis
with sex*diagnosis
. Since diagnosis
has 3 levels itself (nonIBD
, UC
, CD
), this will produce five terms (name
column) for each feature:
-
diagnosisCD
is the difference between female CD and female non-IBD (female non-IBD is the baseline) -
diagnosisUC
is the difference between female UC and female non-IBD -
sexMale
is the difference between male non-IBD and female non-IBD -
sexMale:diagnosisCD
(the interaction ofmale
anddiagnosisCD
) is the difference-in-differences between CD versus non-IBD in males versus females, a measure of whether CD affects males differently from females. That is, given the difference between non-IBD and CD in males and the difference between non-IBD and CD in females, thesexMale:diagnosisCD
coefficient is the difference between the two differences. -
sexMale:diagnosisUC
(the interaction ofmale
anddiagnosisUC
) is the difference-in-differences between UC versus non-IBD in males versus females, a measure of whether UC affects males differently from females.
Another feature of MaAsLin 3 is the ability to test for level-versus-level differences in ordered predictors with contrast tests. Ordered predictors are categorical variables with a natural ordering such as cancer stage, consumption frequency of a dietary factor, or dosage group. Here, we assess how microbial abundances and prevalences are associated with eating red meat by including ordered(red_meat)
in the formula. This will perform a contrast test for whether there is a difference between each pair of subsequent levels (e.g., "Yesterday, 3 or more times" versus "Yesterday, 1 to 2 times") rather than whether there are differences between the levels and the baseline (e.g., "Yesterday, 3 or more times" versus "Not in the last 7 days"). The coefficient, standard error, and p-value all correspond to the difference between the level in value
and the previous level. Ordered predictors should only be included as fixed effects (i.e., no ordered predictors as random effects etc.).
# Put the red meat consumption responses in order
metadata <- metadata %>%
mutate(red_meat = ifelse(red_meat == 'No, I did not consume these products in the last 7 days',
'Not in the last 7 days',
red_meat) %>%
factor(levels = c('Not in the last 7 days',
'Within the past 4 to 7 days',
'Within the past 2 to 3 days',
'Yesterday, 1 to 2 times',
'Yesterday, 3 or more times'))
)
# Create the model with only non-IBD subjects
param_list <- list(input_data = taxa_table,
input_metadata = metadata[metadata$diagnosis == 'nonIBD',],
output = 'ordered_outputs',
formula = '~ ordered(red_meat) + antibiotics + age + reads',
plot_summary_plot = T,
plot_associations = T,
heatmap_vars = c('red_meat Within the past 4 to 7 days',
'red_meat Within the past 2 to 3 days',
'red_meat Yesterday, 1 to 2 times',
'red_meat Yesterday, 3 or more times'),
max_pngs = 30)
fit_out <- maaslin3(param_list)
If we look at the resulting heatmap, we can identify the Alistipes shahii prevalence association as potentially interesting since it increases in prevalence with more meat consumption in all but one level-versus-level comparison.
feature | metadata | value | name | coef | stderr | pval_individual | qval_individual | model | N | N.not.zero | pval_joint | qval_joint |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Alistipes_shahii | red_meat | Within the past 4 to 7 days | red_meatWithin the past 4 to 7 days | 1.47 | 0.353 | 0.0000316 | 0.000701 | logistic | 405 | 266 | 0.0000631 | 0.00113 |
Alistipes_shahii | red_meat | Within the past 2 to 3 days | red_meatWithin the past 2 to 3 days | -0.67 | 0.316 | 0.0341 | 0.164 | logistic | 405 | 266 | 0.067 | 0.22 |
Alistipes_shahii | red_meat | Yesterday, 1 to 2 times | red_meatYesterday, 1 to 2 times | 1.07 | 0.291 | 0.000241 | 0.00411 | logistic | 405 | 266 | 0.000483 | 0.00649 |
Alistipes_shahii | red_meat | Yesterday, 3 or more times | red_meatYesterday, 3 or more times | 3.29 | 3.83 | 0.391 | 0.652 | logistic | 405 | 266 | 0.00677 | 0.0484 |
The last feature of MaAsLin 3 highlighted here is the ability to test for group-wise differences in categorical predictors using an ANOVA or ANOVA-like procedure. Group-wise predictors are categorical variables with or without an ordering such as race, country, or consumption frequency of a dietary factor. When performing a group-wise difference test, we are not comparing any two particular levels; rather, we are investigating whether abundances and prevalences are the same for all levels. Therefore, no coefficient is returned, only a p-value. Here, we test whether there are differences in microbial abundances and prevalences across people with different levels of red meat consumption by including group(red_meat)
. Group-wise predictors should only be included as fixed effects (i.e., no group-wise predictors as random effects etc.).
param_list <- list(input_data = taxa_table,
input_metadata = metadata[metadata$diagnosis == 'nonIBD',],
output = 'group_outputs',
formula = '~ group(red_meat) + antibiotics + age + reads',
plot_summary_plot = T,
plot_associations = T,
heatmap_vars = c('red_meat Within the past 4 to 7 days',
'red_meat Within the past 2 to 3 days',
'red_meat Yesterday, 1 to 2 times',
'red_meat Yesterday, 3 or more times'),
max_pngs = 200)
fit_out <- maaslin3(param_list)
If we look at the same Alistipes shahii association from before, we see that no coefficient or standard error is returned for a group predictor, but the q-value is very low. This corroborates the observation earlier that there are differences in Alistipes shahii prevalence with red meat consumption.
feature | metadata | value | name | coef | stderr | pval_individual | qval_individual | model | N | N.not.zero | pval_joint | qval_joint |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Alistipes_shahii | red_meat | NA | red_meat | NA | NA | 0.00000000586 | 0.000000342 | logistic | 405 | 266 | 0.0000000117 | 0.00000057 |
MaAsLin 3 can also be run with a command line interface. For example, the first HMP2 analysis can be performed with:
./R/maaslin3.R inst/extdata/HMP2_taxonomy.tsv inst/extdata/HMP2_metadata.tsv command_line_output --formula='~ diagnosis + dysbiosis_state + antibiotics + age + reads' --reference='diagnosis,nonIBD;dysbiosis_state,none;antibiotics,No'
- Make sure to provide the full path to the MaAsLin3 executable (i.e.
./R/maaslin3.R
). - In the demo command:
-
inst/extdata/HMP2_taxonomy.tsv
is the path to your data (or features) file -
inst/extdata/HMP2_metadata.tsv
is the path to your metadata file -
command_line_output
is the path to the folder to write the output
-
- HUMAnN 2.0
- HUMAnN 3.0
- MetaPhlAn 2.0
- MetaPhlAn 3.0
- MetaPhlAn 4.0
- MetaPhlAn 4.1
- PhyloPhlAn 3
- PICRUSt 2.0
- ShortBRED
- PPANINI
- StrainPhlAn 3.0
- StrainPhlAn 4.0
- MelonnPan
- WAAFLE
- MetaWIBELE
- MACARRoN
- FUGAsseM
- HAllA
- HAllA Legacy
- ARepA
- CCREPE
- LEfSe
- MaAsLin 2.0
- MaAsLin 3.0
- MMUPHin
- microPITA
- SparseDOSSA
- SparseDOSSA2
- BAnOCC
- anpan
- MTXmodel
- MTX model 3
- PARATHAA