fit: Fit model and test for differential expression

In bhagwataditya/importomics: Unified statistal Modeling of Omics Data

Fit model and test for differential expression

Description

Fit model and test for differential expression

Usage

fit(
  object,
  formula = default_formula(object),
  engine = "limma",
  drop = varlevels_dont_clash(object, all.vars(formula)),
  codingfun = contr.treatment,
  design = create_design(object, formula = formula, drop = drop, codingfun = codingfun),
  contrasts = NULL,
  coefs = if (is.null(contrasts)) colnames(design) else NULL,
  block = NULL,
  weightvar = if ("weights" %in% assayNames(object)) "weights" else NULL,
  statvars = c("effect", "p", "se", "t")[1:2],
  sep = FITSEP,
  suffix = paste0(sep, "limma"),
  verbose = TRUE,
  plot = FALSE
)

fit_limma(
  object,
  formula = default_formula(object),
  drop = varlevels_dont_clash(object, all.vars(formula)),
  codingfun = contr.treatment,
  design = create_design(object, formula = formula, drop = drop, codingfun = codingfun),
  contrasts = NULL,
  coefs = if (is.null(contrasts)) colnames(design) else NULL,
  block = NULL,
  weightvar = if ("weights" %in% assayNames(object)) "weights" else NULL,
  statvars = c("effect", "p"),
  sep = FITSEP,
  suffix = paste0(sep, "limma"),
  verbose = TRUE,
  plot = FALSE
)

.fit_limma(
  object,
  formula = default_formula(object),
  drop = varlevels_dont_clash(object, all.vars(formula)),
  codingfun = "treatment",
  design = create_design(object, formula = formula, drop = drop, codingfun = codingfun),
  contrasts = NULL,
  coefs = if (is.null(contrasts)) colnames(design) else NULL,
  block = NULL,
  weightvar = if ("weights" %in% assayNames(object)) "weights" else NULL,
  statvars = c("effect", "p", "se", "t")[1:2],
  sep = FITSEP,
  suffix = paste0(sep, "limma"),
  verbose = TRUE,
  plot = FALSE
)

fit_wilcoxon(
  object,
  formula = default_formula(object),
  drop = NULL,
  codingfun = contr.treatment.explicit,
  contrasts = NULL,
  coefs = NULL,
  block = NULL,
  weightvar = NULL,
  statvars = c("effect", "p"),
  sep = FITSEP,
  verbose = TRUE,
  plot = FALSE
)

Arguments

object	SummarizedExperiment
formula	modeling formula
engine	'limma', 'lm', 'lme', 'lmer', or 'wilcoxon'
drop	TRUE or FALSE
codingfun	factor coding function contr.treatment: intercept = y0, coefi = yi - y0 contr.treatment.explicit: intercept = y0, coefi = yi - y0 code_control: intercept = ymean, coefi = yi - y0 contr.diff: intercept = y0, coefi = yi - y(i-1) code_diff: intercept = ymean, coefi = yi - y(i-1) code_diff_forward: intercept = ymean, coefi = yi - y(i+) code_deviation: intercept = ymean, coefi = yi - ymean (drop last) code_deviation_first: intercept = ymean, coefi = yi - ymean (drop first) code_helmert: intercept = ymean, coefi = yi - mean(y0:(yi-1)) code_helmert_forward: intercept = ymean, coefi = yi - mean(y(i+1):yp)
design	design matrix
contrasts	NULL or character vector: coefficient contrasts to test
coefs	NULL or character vector: model coefs to test
block	block svar (or NULL)
weightvar	NULL or name of weight matrix in assays(object)
statvars	character vector: subset of c('effect', 'p', 'fdr', 't', 'se')
sep	string: pvar separator ("~" in "p~t2~limma")
suffix	string: pvar suffix ("limma" in "p~t2~limma")
verbose	whether to msg
plot	whether to plot

Value

Updated SummarizedExperiment

Examples

# Read
    file <- system.file('extdata/atkin.metabolon.xlsx', package = 'autonomics')
    object <- read_metabolon(file)
    
# Standard
    fdt(object) %<>% extract(, 'feature_id')
    object %<>% fit_lm(        ~ subgroup)                     #     statistics default
    object %<>% fit_limma(     ~ subgroup)                     # bioinformatics default
    summarize_fit(fdt(object))
    
# Blocked
    fdt(object) %<>% extract(, 'feature_id')
    object %<>% fit_limma(     ~ subgroup, block = 'Subject')  #        simple random effects
    object %<>% fit_lme(       ~ subgroup, block = 'Subject')  #      powerful random effects
    object %<>% fit_lmer(      ~ subgroup, block = 'Subject')  # more powerful random effects
    summarize_fit(fdt(object))
    
# Alternative coding: e.g. grand mean intercept
    fdt(object) %<>% extract(, 'feature_id')
    object %<>% fit_limma(     ~ subgroup, block = 'Subject', codingfun = code_control)
    object %<>% fit_lme(       ~ subgroup, block = 'Subject', codingfun = code_control)
    object %<>% fit_lmer(      ~ subgroup, block = 'Subject', codingfun = code_control)
    summarize_fit(fdt(object))
    
# Posthoc contrasts (only limma!)
    fdt(object) %<>% extract(, 'feature_id')
    object %<>% fit_limma(     ~ subgroup, block = 'Subject', codingfun = code_control, coefs = 't1-t0')
    object %<>% fit_limma( ~ 0 + subgroup, block = 'Subject', contrasts = 't1-t0')
        # flexible, but only approximate
        # stat.ethz.ch/pipermail/bioconductor/2014-February/057682.html
        
# Non-parametric: wilcoxon
    fdt(object) %<>% extract(, 'feature_id')
    object %<>% fit_wilcoxon( ~ subgroup)                    # unpaired
    object %<>% fit_wilcoxon( ~ subgroup, block = 'Subject') # paired
    
# Custom separator
    fdt(object) %<>% extract(, 'feature_id')
    fdt( fit_lm(      object, sep = '.'))
    fdt( fit_limma(   object, block = 'Subject', sep = '.') )
    fdt( fit_lme(     object, block = 'Subject', sep = '.') )
    fdt( fit_lmer(    object, block = 'Subject', sep = '.') )
    fdt( fit_wilcoxon(object, block = 'Subject', sep = '.') )
    fdt( fit_wilcoxon(object, sep = '.') )

bhagwataditya/importomics documentation built on May 1, 2024, 2:01 a.m.