Skip to contents

Formula-based Linear Models

gly_linear_model.Rd

Fit a linear model to every variable in a glycomics or glycoproteomics experiment using limma's empirical Bayes moderation.

Usage

gly_linear_model(
  exp,
  formula,
  contrasts = NULL,
  p_adj_method = "BH",
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

exp

A glyexp::GlycomicSE() or glyexp::GlycoproteomicSE() object, or another SummarizedExperiment containing expression data and sample information.

formula

A one-sided formula describing the model in terms of columns in the sample information, for example ~ treatment * time + batch.

contrasts

An optional uniquely named character vector of limma-style contrast expressions. Use backticks around non-syntactic coefficient names, such as `treatmentB:timeT2`.

p_adj_method

A character string specifying the method used to adjust p-values within each coefficient or contrast. See stats::p.adjust() for available methods. If NULL, no adjustment is performed.

add_info

A logical value. If TRUE (default), variable information from the experiment is added to the tidy result.

...

Additional arguments passed to limma::lmFit().

log_transform

A single non-missing logical value. If TRUE (default), transform expression values using log2(x + 1e-6) before analysis. If FALSE, use the supplied values unchanged, for example for already log-transformed data.

Value

A list with classes glystats_linear_model_res and glystats_res containing:

  • tidy_result: A tibble with variable, term, estimate, AveExpr, t, p_val, optional p_adj, and b columns.

  • raw_result: A list containing the moderated limma fit, design matrix, contrast matrix, and coefficient-name mapping.

  • meta_data: Metadata copied from exp.

Details

gly_linear_model() and gly_limma() both use limma under the hood, but serve complementary use cases. gly_linear_model() provides a general formula interface for specifying a wide range of analysis designs and custom contrasts. gly_limma() is a dedicated wrapper for common differential expression analysis tasks, with convenient group, covariate, subject, and pairwise-comparison arguments.

By default, expression values are transformed using log2(x + 1e-6) before modeling. Set log_transform = FALSE to use the supplied values unchanged. Coefficient estimates and average expression are reported on the selected analysis scale. When contrasts is NULL, all non-intercept coefficients are tested. When contrasts are supplied, only the named contrasts are returned.

Character predictors are handled as factors by stats::model.matrix(). Factor reference levels therefore determine the meaning of coefficients. Set factor levels in the sample information before calling this function when a specific reference level is required.

P-values are adjusted independently across variables for each coefficient or contrast.

See also

gly_limma(), limma::lmFit(), limma::eBayes()

Examples

exp <- glyexp::real_experiment
complete <- rowSums(is.finite(SummarizedExperiment::assay(exp))) == ncol(exp)
exp <- exp[which(complete)[seq_len(10)], ]
SummarizedExperiment::colData(exp)$time <- factor(
  rep(c("early", "late"), length.out = ncol(exp))
)

model_res <- gly_linear_model(exp, ~ group * time)

contrast_res <- gly_linear_model(
  exp,
  ~ group * time,
  contrasts = c(
    C_late_vs_H_late = "groupC + `groupC:timelate`"
  )
)