Hi all,

I'm doing WGCNA downstream to DESeq2 and would like to correct for the effect of covariates, including continuous ones (SVs from using SVA package, pH, Age, PMI) and a discrete covariate, namely sex, before moving on with transforming counts using VST. Importantly, I have no batches to correct for.

Below is the code I used when I didn't account for covariates.

dds <- DESeqDataSetFromMatrix(countData = rawcounts, colData = sample.info, design = ~ SV1 + SV2 + pH + Age + Sex + PMI + Group)
transformed.counts <- varianceStabilizingTransformation(dds, blind = TRUE) # try blind = FALSE if wanting to include design formula
transformed.counts <- assay(transformed.counts)

What is the best way to correct for the mentioned covariates? I have read this post proposing use of empiricalBayesLM, but have come across other people using limma::removeBatchEffect. Which one is better for my purpose, if any? admittedly I find the code from the limma::removeBatchEffect post easier to follow and adapt to my purpose. Could anyone help me out with the code needed to do the job?

Many thanks!