MaAsLin 3 Tutorial

Running MaAsLin 3 on HMP2 data

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', base 2) 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 (see normalization and transform 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 (see augment 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 (see max_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 parameters.
• 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 need not be filtered out since the prevalence modeling in MaAsLin 3 already accounts for high proportions of zeros. However, filtering low prevalence features might improve power.

set.seed(1)
fit_out <- maaslin3(input_data = taxa_table,
                    input_metadata = metadata,
                    output = 'hmp2_output',
                    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 = 10,
                    cores = 1)

By default, many lines of information are printed while the models fit, but the verbosity can be changed with the verbosity argument (e.g., verbosity = 'WARN').

MaAsLin 3 can also take as input a SummarizedExperiment or TreeSummarizedExperiment object from BioConductor:

se <- SummarizedExperiment(
    assays = list(taxa_table = t(taxa_table)),
    colData = metadata
)

fit_out <- maaslin3(input_data = se,
                    output = 'hmp2_output',
                    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 = 10,
                    cores = 1)

Median comparisons

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 is TRUE by default and should be used to make inference on the absolute scale when using relative abundance data. When median_comparison_abundance is TRUE, only the p-values and q-values change; the coefficients returned are still the relative abundance coefficients unless subtract_median is TRUE in which case the median will be subtracted from the coefficients.
• median_comparison_abundance should be FALSE 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 with unscaled_abundance.
• median_comparison_prevalence is FALSE by default.

Significant associations

The main output from MaAsLin 3 is the list of significant associations in significant_results.tsv. This lists all associations with joint or individual q-values that pass MaAsLin 3’s significance threshold, ordered by increasing individual q-value. The format is:

• feature and metadata are the feature and metadata names.
• value and name are the value of the metadata and variable name from the model.
• coef and stderr 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^coef fold change in the relative abundance of the feature. In prevalence models, a one-unit change in the metadatum variable corresponds to a coef change in the log-odds of a feature being present.
• pval_individual and qval_individual are the p-value and false discovery rate (FDR) corrected q-value of the individual association. FDR correction is performed over all p-values without errors in the abundance and prevalence modeling together.
• pval_joint and qval_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.
• model specifies whether the association is abundance or prevalence.
• N and N_not_zero are the total number of data points and the total number of non-zero data points for the feature.

Full output file structure

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).

1. 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_linear.rds and models_logistic.rds
     • These files contain a list with every model fit object (linear for linear models, logistic for logistic models).
     • It will only be generated if save_models is set to TRUE.
   • residuals_linear.rds and residuals_logstic.rds
     • These files contain a data frame with residuals for each feature.
   • fitted_linear.rds and fitted_logistic.rds
     • These files contain a data frame with fitted values for each feature.
   • ranef_linear.rds and ranef_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.
2. 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 below max_significance divided by 10. (Make sure to set include=T to view the plots when knitting the Rmd.)
   • association_plots/[metadatum]/[association]/ [metadatum]_[feature]_[association].png
     • 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.
     • To see the association plots, set max_pngs = 250 above and set eval=TRUE and include=TRUE in this chunk

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.

Diagnostics

There are a few common issues to check for in the results:

1. When warnings or errors are thrown during the fitting process, they are recorded in the error column of all_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.
   • With warn_prevalence=TRUE, a note will be produced in the error column indicating when prevalence associations are likely induced by abundance associations. These should be checked visually with the diagnosic visuals.
2. 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.
3. The q-values are FDR corrected over all abundance and prevalence relationships, so it may be preferable to FDR correct just the p-values from the variables of interest. This can reduce false positives when there are many significant but uninteresting associations (e.g., many read depth associations).
4. 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.

Replotting

Once maaslin3 has been run once, maaslin_plot_results_from_output can be 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 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
heatmap_vars = c('Dysbiosis UC', 'Diagnosis UC',
                'Abx Used', 'Age', 'Read depth')
coef_plot_vars = c('Dysbiosis CD', 'Diagnosis CD')

# This section is necessary for updating the association plots
taxa_table_copy <- taxa_table
colnames(taxa_table_copy) <- gsub('_', ' ', colnames(taxa_table_copy)) %>%
    gsub(pattern = 'sp ', replacement = 'sp. ')

# Rename the features in the norm transformed data file
data_transformed <-
    read.csv('hmp2_output/features/data_transformed.tsv', sep='\t')
colnames(data_transformed) <-
    gsub('_', ' ', colnames(data_transformed)) %>%
    gsub(pattern = 'sp ', replacement = 'sp. ')
write.table(data_transformed,
            'hmp2_output/features/data_transformed.tsv',
            sep='\t', row.names = FALSE)

# Rename the metadata like in the outputs table
metadata_copy <- metadata
colnames(metadata_copy) <-
    case_when(colnames(metadata_copy) == 'age' ~ 'Age',
            colnames(metadata_copy) == 'antibiotics' ~ 'Abx',
            colnames(metadata_copy) == 'diagnosis' ~ 'Diagnosis',
            colnames(metadata_copy) == 'dysbiosis_state' ~ 'Dysbiosis',
            colnames(metadata_copy) == 'reads' ~ 'Read depth',
            TRUE ~ colnames(metadata_copy))
metadata_copy <- metadata_copy %>%
    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(
    output = 'hmp2_output',
    metadata = metadata_copy,
    normalization = "TSS",
    transform = "LOG",
    median_comparison_abundance = TRUE,
    median_comparison_prevalence = FALSE,
    max_significance = 0.1,
    max_pngs = 20)

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.

Spike-in

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]

##            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

spike_in_metadata[1:5,]

##         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

spike_in_unscaled[1:5, , drop=FALSE]

##         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.

fit_out <- maaslin3(
    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)

The absolute abundance coefficients can be inspected:

rownames(fit_out$fit_data_abundance$results) <- NULL
head(fit_out$fit_data_abundance$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")

Total abundance scaling

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]

##          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

total_scaling_metadata[1:5,]

##         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_scaling_unscaled[1:5, , drop=FALSE]

##             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.

fit_out <- maaslin3(
    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)

The absolute abundance coefficients can be inspected:

rownames(fit_out$fit_data_abundance$results) <- NULL
head(fit_out$fit_data_abundance$results, n = 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")

Computationally estimated absolute abundance

When using a purely computational approach to estimate absolute abundance (e.g., Nishijima et al. Cell 2024 with the function MLP), the abundance table can be estimated and then run with the normalization = 'NONE' to avoid scaling the results.