Hi everyone,

I am analyzing Bulk RNA-Seq data, I generated the counts with HTSeqcount and I analyzed them with DESeq2.

Then, I want to get the normalized counts corrected for the variables in my design matrix (Like age, sex, RIN, ...) to use it as an input for WGCNA.

But I realized when I am using variance normalized counts (vst-counts) from DESeq2 they are not corrected for variables in my design matrix.

I was wondering which normalized counts I can use that can be corrected for variables in my design matrix, therefore it can be used as an input for WGCNA?

I also tried to use "limma::removeBatchEffect" but still, it didn't really correct my count for those variables.

I would be grateful if anyone could provide/recommend any solutions, please!

Here is the code that I used in my DESeq2:

dds <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable,
                                  directory=folder,
                                  design=~Plate+RIN+Sex+Age+condition+PC3+PC2+PC1)

dds <- estimateSizeFactors(dds)
keep <- rowSums( counts(dds) >= 10 ) >= 20 #filters out genes 
dds <- dds[keep,]

colData(dds)$condition <- relevel(colData(dds)$condition, ref = "Control")

dds<- DESeq(dds)

vsd <- vst(dds, blind = FALSE)

I want to get the normalized counts adjusted for sex, RIN, plate, age, PC1, PC2, PC3

Thank you in advance!