I'm aware of that limma functions were originally designed for microarray and RNA-Seq data, but I'm using them for modeling associations between metabolic measures (quantified with nuclear magnetic resonance (NMR) platform) and different traits (e.g. body mass index and overweight) here. I prefer using lmFit() and eBayes() over standard regression functions like lm() because, as far as I have understood, these limma functions account for shared structure of in high-dimensional data, optimize inference across multiple tests and moderate the variance estimates with the empirical Bayes approach.

I'm running the following chunk of code inside a loop over different traits some of which are binary and others are continuous. The aim is to test associations between each metabolic measure and trait.

# ph.dat is a data frame with individuals in rows, a single column for individual IDs, and additional columns for traits
# mtb.dat is a matrix with metabolic measures in rows and individuals (the same that are in ph.dat) in columns
# gr.labs is a named list where each element is a vector including labels (e.g. 'cases', 'controls') for a specific binary trait

stopifnot(all(ph.dat$id == colnames(mtb.dat)))

if(trait %in% bin.traits) { # if trait is binary
  gr.labs0 = gr.labs[[trait]]
  groupyes <- which(ph.dat[,trait]==gr.labs0[1])
  groupno <- which(ph.dat[,trait]==gr.labs0[2])

  design <- model.matrix(~0 + ph.dat[,trait])
  colnames(design) <- gsub("ph.dat\\[, trait\\]", "", colnames(design))
  design <- design[,c(gr.labs0[2], gr.labs0[1])]
  colnames(design) <- c("groupno", "groupyes")
  cont.matrix <- makeContrasts(yesvsno = (groupyes - groupno), levels = design)
  fit <- lmFit(mtb.dat, design)
  fit2 <- contrasts.fit(fit, cont.matrix)
} else { # if trait is not binary, it is assumed to be continuous
  design <- model.matrix(~ ph.dat[,trait])
  fit2 <- lmFit(mtb.dat, design)
}

fit.e <- eBayes(fit2)
topDM <- topTable(fit.e, number=nrow(mtb.dat))

The code runs without errors, but I want to ensure that my approach is conceptually sound. In particular, I would like to confirm the following points:

1. I tried to see what happens inside lmFit() and it seems that lm.series() is applied regardless of whether the trait is continuous or binary. The association is then modeled so that the trait is the predictor and 'omical' variable (here metabolic measure) is the response. Does lm.series() make some assumptions on the type or structure of the omics data? Is there some specific in lm.series() that would make it unsuitable for modeling associations with metabolic measures? What are its key differences compared to e.g. lm()? Is lm.series() just a wrapper around running e.g. lm() repeatedly for each feature (row)?
2. Do lmFit(), eBayes() or topTable() make any assumptions about whether the trait is continuous or categorical?
3. What is the purpose of the trend argument in eBayes(), and under what circumstances is it recommended to set it to TRUE versus leaving it at the default FALSE? According to the limma user guide (http://bioconductor.jp/packages/3.9/bioc/vignettes/limma/inst/doc/usersguide.pdf) setting trend = TRUE may increase statistical power but for example in the manual https://bioconductor.org/packages/release/bioc/manuals/limma/man/limma.pdf this argument is left at its default FALSE.

The assumptions made by limma are explained in the documentation and in the relevant publications:

• https://doi.org/10.1214/16-aoas920
• https://doi.org/10.1093/nar/gkv007
• https://gksmyth.github.io/pubs/ebayes.pdf

limma works with any number of binary, categorical or continuous predictors. It works exactly the same in each case.

Regarding the trend argument, try it and see. If there is clear trend then it will strengthen the analysis.