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, input_metadata, 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[] |> 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:
That worked, and this is exactly what I was looking for. Thanks so much for your help!