I have recently tried to implement the Samon-Tximport-DESeq2 pipeline with my latest RNA-Seq dataset however I have a problem with some of the downstream analysis and potential compatibility with other programs. I have found that the dataset has quite large batch effects and I am hoping to regress these out prior to downstream analysis (mainly WGCNA). My general workflow at present has been to use tximport to get gene-level count estimates and then to run the standard DESeq2 pipeline from there.

I was then planning to normalize these counts (vst) and run a simple regression on the transformed object to remove the effects of batch. The residuals I would then use for WGCNA. My concern however is that it is mentioned throughout the documentation that one should not normally use these transformed objects for any purpose involving DE. I was hoping to get some advice on whether this would be appropriate and if not, what other methods may be more suitable. I have a general population cohort and I am looking at more covariates of interest than disease status (age, gender, family history etc etc).

Many thanks,

Matthew

You normally don't use residuals from a model fit as input to another model fit, because you are then cheating by ignoring the fact that you have lost degrees of freedom (in the original model fit) that you are not accounting for in the second model fit. This isn't an issue for things like WGCNA, which aren't inferential techniques. In other words, you won't be using the residuals for differential expression, but instead will be using for WGCNA, which doesn't do that sort of thing.