This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Batch correction in DESeq2---Normalized Counts

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!

rna-seq
• 4,486 views • link
updated by ATpoint ATpoint • written by mtsompana • mtsompana

If you are looking for gene level comparison consider only the normalized values don't take rlog

• 0 views • link
written by DareDevil DareDevil

1 answer

If you want to start with the normalized counts without the extended rlog procedure on top of it (rlog are normalized counts plus variance-stabilized and fold changes being moderated/shrunken, see manual) then simply use:

normbatch <- limma::removeBatchEffect(log2(counts(dds, normalized=TRUE)+1), dds$Batch)
• 0 views • link
written by ATpoint ATpoint

Thank you so much. Do you recommend using the rlog values instead?

• 0 views • link
written by mtsompana • mtsompana

I think I figured it out. Could you please confirm this is correct. Thank you!!

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ Batch + Condition3)

dds <- dds[ rowSums(counts(dds)) > 10, ]

dds <- DESeq(dds)

rld<-rlog(dds,blind = FALSE)

assay(rld) <- limma::removeBatchEffect(assay(rld), rld$Batch)

counts_afterbatchcorrection <- assay(rld)

write.csv(cts.norm, file = "normalizedCountsafterbatchcorrection.csv")

• 0 views • link
written by mtsompana • mtsompana

Yeah, this is the same as in your original question, is it? Looks ok to me. rlog has a couple of nice features and it is recommended by the DESeq2 author for downstream applications. I use it a lot.

• 0 views • link
written by ATpoint ATpoint

thank you for your time and help! yes this is for the same question.

• 0 views • link
written by mtsompana • mtsompana

Log in to answer this question.

More posts like this