Hello,

I have a RNA-seq dataset in which I am adjusting for subject age (design=~age+condition).

I am wondering if DESeq2 generates a table of normalized reads after adjusting for the covariate, "age".  If so, how can I extract it?

I ask because I would like to extract the normalized reads (adjusted for the covariate) that I can use to create a gene co-expression network.
 

Similar to g.atla's answer, you can perform one of DESeq2's transformations (see vignette) and then, to adjust for a covariate, you can use limma's removeBatchEffect() function on the transformed values. Note that the transformed values will be on the log2 scale.

# you can use VST, rlog, or normTransform() function in DESeq2
# to create a DESeqTransform object. then:
mat <- assay(vsd) 
library(limma)
mat2 <- removeBatchEffect(mat, vsd$age.binned)

What's this age.binned I use in the last line? Read the DESeq2 FAQ about using age as a continuous covariate. In short, I do not recommend adjusting for age as a continuous covariate, and offer an alternative by creating a new variable which bins age into, say, 5 groups using the cut() function.

vsd$age.binned <- cut(vsd$age, 5)