Intercept value in maaslin2

Hi! I’m trying to obtain the value of the intercept for each of the features in my MaAsLin2 results. Is it possible to generate a column for this in the output? I have fitted a few models with interaction terms where this would be especially helpful for following up with additional analyses and plots. Thanks!

If you set Maaslin2(..., save_models = TRUE), among the outputs you’ll get an rds file containing the model fits. You can extract the intercepts from there. Using the example in the package:

fit_data <- Maaslin2(input_data, 
                     output = '~/tmp', 
                     transform = "AST",
                     fixed_effects = c('diagnosis', 'dysbiosisnonIBD','dysbiosisUC','dysbiosisCD', 'antibiotics', 'age'),
                     random_effects = c('site', 'subject'),
                     normalization = 'NONE',
                     reference = 'diagnosis,nonIBD',
                     standardize = FALSE, 
                     save_models = TRUE) 

models <- readRDS("~/tmp/fits/models.rds")

models[[1]] |> summary() |> coefficients() 

#                        Estimate   Std. Error        df    t value    Pr(>|t|)
#(Intercept)         2.655611e-02 0.0092456815  121.0273  2.8722722 0.004813519
#diagnosisCD        -1.812594e-02 0.0088601782  125.6027 -2.0457758 0.042866199
#diagnosisUC         4.091809e-03 0.0096240661  123.9921  0.4251643 0.671453859
#dysbiosisnonIBDYes -6.719608e-03 0.0071908428 1551.8470 -0.9344674 0.350208208
#dysbiosisUCYes     -1.056876e-02 0.0054566951 1498.1532 -1.9368420 0.052952474
#dysbiosisCDYes      6.020841e-04 0.0041147376 1567.5584  0.1463238 0.883684589
#antibioticsYes     -6.128646e-03 0.0033373554 1526.8368 -1.8363779 0.066496111
#age                 4.813489e-05 0.0001984042  122.9060  0.2426102 0.808711741

Thanks for your quick response! I just tried that and unfortunately received this error:

Error in Maaslin2(input_data = asv.relabun, input_metadata = input_metadata,  : 
  unused argument (save_models = TRUE)

In case it’s helpful, this is what I am running. It normally works fine without the “save_models = TRUE” parameter:

mod3 <- Maaslin2(
  input_data = asv.relabun, input_metadata = input_metadata, output = "output_file",
  fixed_effects = c("timepoint_type","GroupHL_01","GroupHL_int"),
  random_effects = "participant_id",
  reference = "timepoint_type,Pre",
  normalization = "NONE",
  save_models = TRUE)

Oops, I didn’t realize that option hadn’t made it to the BioC version yet. You can get it by installing the development version from github: devtools::install_github("biobakery/Maaslin2")

That worked, and this is exactly what I was looking for. Thanks so much for your help!

1 Like