Hi everyone, I'm performing analysis of some RNAseq samples, and currently trying to cope with batch effect.plotting PCA of vsd transformed data, I can clearly see two batches which are differ fromt the others. 

    plotPCA(vsd, intgroup = c('batch'),ntop= 34085)
![PCA][1]

If I'm not mistaken, then for DE expression analysis I could use design formula of dds function, to reduce batch effect:

    dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design =~ batch + type)
    dds <- DESeq(dds)
    vsd <- varianceStabilizingTransformation(dds, blind = TRUE)


But I also would like to visualise results of batch correction. For PCA plot I could simply follow the DESeq2 manual and just perform limma batch effect removal:

    mat <- assay(vsd)
    mat <- limma::removeBatchEffect(mat, vsd$batch)
    assay(vsd) <- mat
    plotPCA(vsd, intgroup = c('batch'),ntop= 34085)

Which gives me a nice results:
![batch corrected PCA][2]

But now I'd like to visualise expression of individual genes of interest. I could use vsd transforemd data, but these numbers mask the actual count numbers, therefore laking info on actual level of gene expression, which I'd like to preserve.
![expression of gene X fo given conditions, vsd transformed data, limma correction][3]

So what would be the better solution here? Could I use limma removeBatchEffect() on normalised count table, maybe? 

Thanks! 

  [1]: https://i.ibb.co/bzspVz1/Intitial-pca.png
  [2]: https://i.ibb.co/RQhzx2x/batch-corrected-PCA.png
  [3]: https://i.ibb.co/Cs6jD0j/geneX.png
`vst` transforms to approximately log2-scale, that is simply a data transformation to allow a wide range of counts to fit on a narrow scale rather than having counts spread between 0 and numbers in the hundreds-of- thousands.