I performed DESEQ on my data implementing the design = ~ condition + batch. I got the result below.

I then corrected the batch effect with Limma and was able to get a better PCA result:

Do I use the Limma batch corrected values to perform my downstream analysis or do I simply proceed with the original DESEQ data? I am reading about it in the forums and it says it is not recommended, but I do not want to proceed to downstream analysis with this batch effect present.

I truly will appreciate any insight. Thank you!

There are two separate paths in this workflow; the one we will see first involves transformations of the counts in order to visually explore sample relationships. In the second part, we will go back to the original raw counts for statistical testing. This is critical because the statistical testing methods rely on original count data (not scaled or transformed) for calculating the precision of measurements.

https://bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#exploratory-analysis-and-visualization

For controlling for batch during statistical testing you would use a multi-factor design with the original counts (not regressed, transformed counts), following the standard workflow.

E.g. here at the top of the vignette:

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#quick-start

or in more depth here:

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs

removeBatchEffect is just for visualization (and its only an approximation -- visualizing batch adjusted transformed counts is not the same as what happens in a multi-factor GLM).

It is possible to visualize the transformed data with batch variation removed, using the removeBatchEffect function from limma.

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#why-after-vst-are-there-still-batches-in-the-pca-plot