I am analyzing bulk RNA-seq data from samples sequenced at different time points (different library preparation and sequencing strategy).
I am using the code below to accommodate for batch effects so I can do clustering:
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ Batch + Condition)
dds <- dds[ rowSums(counts(dds)) > 10, ]
rld<-rlog(dds,blind = FALSE)
assay(rld) <- limma::removeBatchEffect(assay(rld), rld$Batch)
I need to get *normalized gene counts (not rlog counts) after I correct for batch effects. I need this so I can compare levels of expression for a specific gene among the different samples. I am not sure how to do this.*
I know how to do it with the original dds object:
dds.normCts <- estimateSizeFactors(dds)
cts.norm <- counts(dds.normCts,normalized=TRUE)
write.csv(cts.norm, file = "normalizedCounts.csv")
Thank you for your time and help in advance!