Skip to contents

Compute linear models (LM) on microbiome data

Source: R/run_lm_model.R
run_lm_model.Rd

This function takes an animalcules-formatted MultiAssayExperiment and runs an independent linear model for each taxon. The model predicts taxon log CPM abundance as a product of user-specified covariates. This model can be used for general microbiome analyses without repeated measures data. Note, the "broom", and "broom.mixed" packages are required to use the testing functionality of this package; both can be installed via CRAN.

Usage

run_lm_model(
  dat,
  taxon_level = "genus",
  fixed_cov,
  p_adj_method = "fdr",
  plot_out = FALSE,
  plotsave_loc = ".",
  plot_terms = NULL,
  ...
)

Arguments

dat

A MultiAssayExperiment object specially formatted as an animalcules output.

taxon_level

Character string, default is "genus".

fixed_cov

A character vector naming covariates to be tested.

p_adj_method

A character string specifying the correction method. Can be abbreviated. See details. Default is "fdr".

plot_out

Logical indicating whether plots should be output alongside the model results. Default is FALSE.

plotsave_loc

A character string giving the folder path to save plot outputs. This defaults to the current working directory.

plot_terms

Character vector. Which terms should be examined in the plot output? Can overlap with the fixed_cov inputs.

...

Further arguments passed to ggsave for plot creation.

Value

A data.frame of modeling results.

Details

P-values are adjusted for the model coefficients within each taxon. The following methods are permitted: c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")

Examples

dat <- system.file("extdata/MAE.RDS", package = "LegATo") |>
  readRDS() |>
  filter_MAE()
#> The overall range of relative abundance counts between samples is (13218, 3016276) 
#> Number of OTUs that exhibit a relative abundance >3% in at least 5% of the total samples: 17/1690
out <- run_lm_model(dat, fixed_cov = c("timepoint", "HIVStatus"),
                    plot_out = FALSE)
head(out)
#>   Coefficient Coefficient Estimate Standard Error Statistic Unadj p-value
#> 1 (Intercept)         1.583191e+05   38802.912402  4.080084  7.012551e-05
#> 2 (Intercept)         4.151148e+01      19.455765  2.133634  3.435965e-02
#> 3 (Intercept)         1.085153e+01       4.747316  2.285824  2.354463e-02
#> 4 (Intercept)         9.800220e+00       4.584216  2.137818  3.401292e-02
#> 5 (Intercept)         4.237653e+00       2.059365  2.057748  4.119716e-02
#> 6 (Intercept)         5.523260e+00       4.591306  1.202982  2.307169e-01
#>                 Taxon Adj p-value
#> 1               Other 0.001192134
#> 2            Belnapia 0.116725295
#> 3         Actinotalea 0.116725295
#> 4          Leptothrix 0.116725295
#> 5 Limosilactobacillus 0.116725295
#> 6      Staphylococcus 0.392218737