Hello, I had a question about extracting covariate-corrected normalized count values from limma (if possible) for box plots. Essentially, I have some genes which limma labels as differentially expressed when covariates are corrected for in the model, but when I pull the voom-normalized counts to visualize expression differences as boxplots between conditions they do not represent the output of the limma model. I imagine this is because the normalized counts are not corrected for any covariates as they are in the model, which is why the differences in expression are not apparent. I'm not sure it is needed but here is the relevant code:

mm <- model.matrix(~1, data = dat.abund.norm$samples)

dat.abund.norm.voom <- voomWithQualityWeights(
  dat.abund.norm,
  design=mm,
  plot=TRUE)

mm <- model.matrix(~PREDICTEDCLASS +
                     WBC +
                     Epithelial +
                     SiteID +
                     Plate, 
                   data = dat$targets)

y2 <- lmFit(dat, mm)

fit2 <- eBayes(y2, trend = TRUE)

top.table_WGCNA <- topTable(fit2, coef = c(2,3,4), n = "Inf")

You should use removeBatchEffect for that.