MaAsLin 3 User Manual

Diagnostics

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.
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 or prevalence relationships (separately), 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 plots 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.

In R

The following code identifies associations between patient metadata and microbial species in the HMP2 cohort.

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

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

# Fit models
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 = 20,
                    cores = 1)

The maaslin3 wrapper function runs all steps at once, but equivalent results would be produced from running one step at a time:

maaslin_log_arguments(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 = 20,
                    cores = 1)

read_data_list <- maaslin_read_data(taxa_table,
                                    metadata)
read_data_list <- maaslin_reorder_data(read_data_list$data,
                                    read_data_list$metadata)

data <- read_data_list$data
metadata <- read_data_list$metadata

formulas <- maaslin_check_formula(data,
                                metadata,
                                '~ diagnosis + dysbiosis_state +
                                    antibiotics + age + reads')

formula <- formulas$formula
random_effects_formula <- formulas$random_effects_formula

normalized_data = maaslin_normalize(data,
                                    output = 'hmp2_output',
                                    normalization = 'TSS')

filtered_data = maaslin_filter(normalized_data,
                            output = 'hmp2_output')

transformed_data = maaslin_transform(filtered_data,
                                    output = 'hmp2_output',
                                    transform = 'LOG')

standardized_metadata = maaslin_process_metadata(metadata,
                                                formula,
                                                standardize = TRUE)

maaslin_results = maaslin_fit(filtered_data,
                            transformed_data,
                            standardized_metadata,
                            formula,
                            random_effects_formula,
                            data = data,
                            median_comparison_abundance = TRUE,
                            median_comparison_prevalence = FALSE,
                            cores = 1)

maaslin_write_results(output = 'hmp2_output',
                    maaslin_results$fit_data_abundance,
                    maaslin_results$fit_data_prevalence,
                    random_effects_formula,
                    max_significance = 0.1)

# Invisible to avoid dumping plots
invisible(maaslin_plot_results(output = 'hmp2_output',
                            transformed_data,
                            metadata,
                            maaslin_results$fit_data_abundance,
                            maaslin_results$fit_data_prevalence,
                            normalization = "TSS",
                            transform = "LOG",
                            median_comparison_abundance = TRUE,
                            median_comparison_prevalence = FALSE,
                            max_significance = 0.1,
                            max_pngs = 20))

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
colnames(taxa_table) <- gsub('_', ' ', colnames(taxa_table)) %>%
    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
colnames(metadata) <-
    case_when(colnames(metadata) == 'age' ~ 'Age',
            colnames(metadata) == 'antibiotics' ~ 'Abx',
            colnames(metadata) == 'diagnosis' ~ 'Diagnosis',
            colnames(metadata) == 'dysbiosis_state' ~ 'Dysbiosis',
            colnames(metadata) == 'reads' ~ 'Read depth',
            TRUE ~ colnames(metadata))
metadata <- 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(
    output = 'hmp2_output',
    metadata,
    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. (Make sure to set include=T to view the plots when knitting the Rmd.)

In the association plots below, the taxa and metadata have been renamed to be consistent with the results file from earlier. To see the plots, set max_pngs = 250 above and set eval=TRUE in this chunk.

Contrast testing

We can also perform contrast testing in MaAsLin 3 to check whether a linear combination of the coefficients is significantly different from some constant. Specify a contrast matrix with column names matching the coefficient names and rows with the contrasts of interest. A separate contrast test will be performed for each row.

contrast_mat <- matrix(c(1, 1, 0, 0, 0, 0, 1, 1), 
    ncol = 4, nrow = 2, byrow = TRUE)
    
colnames(contrast_mat) <- c("diagnosisUC",
                            "diagnosisCD",
                            "dysbiosis_statedysbiosis_UC",
                            "dysbiosis_statedysbiosis_CD")
                            
rownames(contrast_mat) <- c("diagnosis_test", "dysbiosis_test")

contrast_mat
                            
maaslin_contrast_test(fits = fit_out$fit_data_abundance$fits,
                        contrast_mat = contrast_mat)

Command line

MaAsLin 3 can also be run with a command line interface. For example, the first HMP2 analysis can be performed with the following command (the slashes may need to be removed):

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

Required parameters

• input_data: A data frame of feature abundances or read counts or a filepath to a tab-delimited file with abundances. It should be formatted with features as columns and samples as rows (or the transpose). The column and row names should be the feature names and sample names respectively.
• input_metadata: A data frame of per-sample metadata or a filepath to a tab-delimited file with metadata. It should be formatted with variables as columns and samples as rows (or the transpose). The column and row names should be the variable names and sample names respectively.
• output: The output folder to write results.

Model formula

• formula: A formula in lme4 format. Random effects, interactions, and functions of the metadata can be included (note that these functions will be applied after standardization if standardize = TRUE). Group, ordered, and strata variables can be specified as: group(grouping_variable), ordered(ordered_variable), and strata(strata_variable). The other variable options below will not be considered if a formula is set.
• fixed_effects: A vector of variable names to be included as fixed effects.
  • Fixed effects models are fit with lm (linear) or glm (logistic).
• reference: For a variable with more than two levels supplied with fixed_effects, the factor to use as a reference provided as a string of ‘variable,reference’ semi-colon delimited for multiple variables.
• random_effects: A vector of variable names to be included as random intercepts. 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 and subsequently filtered out. (See strata_effects as well.)
  • Random effects models are fit with lmer (linear) and glmer (logistic), and the significance tests come from lmerTest and glmer respectively.
• group_effects: A factored categorical variable to be included for group testing. An ANOVA-style test will be performed to assess whether any of the variable’s levels are significant, and no coefficients or individual p-values will be returned.
  • Tests are performed with the anova function’s LRT option (logistic fixed and mixed effects), the anova function’s F test (linear fixed effects), or lmerTest::contest (linear mixed effects).
• ordered_effects: A factored categorical variable to be included. Consecutive levels will be tested for significance against each other with contrast tests, and the resulting associations will correspond to effect sizes, standard errors, and significances of each level versus the previous.
  • Contrast tests are performed with multcomp::glht (fixed effects and logistic mixed effects) and lmerTest::contest (linear mixed effects).
• strata_effects: A single grouping variable to be included in matched case-control studies. If a strata variable is included, no random effects can be included. When a strata variable is included, a conditional logistic regression will be run to account for the strata. The abundance model will be run with a random intercept in place of the strata. Strata can include more than two observations per group. Only variables that differ within the groups can be tested. In general, strata effects are not recommended except for advanced users. Fixed or random intercepts are recommended instead.

Feature specific covariates

Particularly for use in metatranscriptomics workflows, a table of feature-specific covariates can be included. A feature’s covariates will be included like a fixed effect metadatum when fitting the model for that feature. The covariate’s name does not need to be included in the formula.

• feature_specific_covariate: A table of feature-specific covariates or a filepath to a tab-delimited file with feature-specific covariates. It should be formatted with features as columns and samples as rows (or the transpose). The row names and column names should be the same as those of the input_data: the column and row names should be the feature names and sample names respectively.
• feature_specific_covariate_name: The name for the feature-specific covariates when fitting the models.
• feature_specific_covariate_record: Whether to keep the feature-specific covariates in the outputs when calculating p-values, writing results, and displaying plots.

Analysis options

• min_abundance (default 0): Features with abundances of at least min_abundance in min_prevalence of the samples will be included for analysis. The threshold is applied before normalization and transformation.
• min_prevalence (default 0): See above.
• zero_threshold (default 0): Abundances less than or equal to zero_threshold will be treated as zeros. This is primarily to be used when the abundance table has likely low-abundance false positives.
• min_variance (default 0): Features with abundance variances less than or equal to min_variance will be dropped. This is primarily used for dropping features that are entirely zero.
• max_significance (default 0.1): The FDR corrected q-value threshold for significance used in selecting which associations to write as significant and to plot.
• normalization (default TSS): The normalization to apply to the features before transformation and analysis. The option TSS (total-sum scaling) is recommended, but CLR (centered log ratio) and NONE can also be used.
• transform (default LOG, base 2): The transformation to apply to the features after normalization and before analysis. The option LOG is recommended, but PLOG (pseudo-log with a pseudo-count of half the dataset minimum non-zero abundance replacing zeros, particularly for metabolomics data) and NONE can also be used.
• correction (default BH): The correction to obtain FDR-corrected q-values from raw p-values. Any valid options for p.adjust can be used.
• standardize (default TRUE): Whether to apply z-scores to continuous metadata variables so they are on the same scale. This is recommended in order to compare coefficients across metadata variables, but note that functions of the metadata specified in the formula will apply after standardization.
• warn_prevalence (default TRUE): Warn when prevalence associations are likely induced by abundance associations. This requires re-fitting the linear models on the TSS log-transformed data. A prevalence coefficient is flagged if its corresponding abundance coefficient is significantly different from 0, of the same sign, and larger in absolute value.
• augment (default TRUE): To avoid linear separability in the logistic regression, at each input data point, add an extra 0 and an extra 1 observation weighted as the number of predictors divided by two times the number of data points. This is almost always recommended to avoid linear separability while having a minor effect on fit coefficients otherwise.
• evaluate_only (default NULL): To fit only the abundance or only the prevalence models, evaluate_only can be set to abundance or prevalence.

Compositionality corrections

Plotting parameters

• plot_summary_plot (default TRUE): Generate a summary plot of significant associations.
• summary_plot_first_n (default 25): Include the top summary_plot_first_n features with significant associations.
• coef_plot_vars: Vector of variable names to be used in the coefficient plot section of the summary plot. Continuous variables should match the metadata column name, and categorical variables should be of the form "[variable] [level]".

By default, the (up to) two metadata variables with the most significant associations will be plotted in the coefficient plot, and the rest will be plotted in the heatmap. Because predicting the output variable names can be tricky, it is recommended to first run maaslin3 without setting coef_plot_vars or heatmap_vars, look at the names of the variables in the summary plot, and then rerun with maaslin_plot_results_from_output after updating coef_plot_vars and heatmap_vars with the desired variables. * heatmap_vars: Vector of variable names to be used in the heatmap section of the summary plot. Continuous variables should match the metadata column name, and categorical variables should be of the form "[variable] [level]". * plot_associations (default TRUE): Whether to generate plots for significant associations. * max_pngs (default 30): The top max_pngs significant associations will be plotted.

Technical parameters

• cores (default 1): How many cores to use when fitting models. (Using multiple cores will likely be faster only for large datasets or complex models.)
• save_models (default FALSE): Whether to return the fit models and save them to an RData file.
• verbosity (default 'FINEST'): The level of verbosity for the logging package.
• save_plots_rds (default FALSE): Whether to save the plots as RDS files.