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)
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:
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.
So what would be the better solution here? Could I use limma removeBatchEffect() on normalised count table, maybe?