Question: Individual genes plots after batch correction with Deseq2
gravatar for Eugene A
20 months ago by
Eugene A90
Eugene A90 wrote:

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: batch corrected PCA

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

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


rna-seq deseq2 R • 1.4k views
ADD COMMENTlink modified 20 months ago • written 20 months ago by Eugene A90

Please see How to add images to a Biostars post to add your images properly. You need the direct link to the image, not the link to the webpage that has the image embedded (which is what you have used here)

ADD REPLYlink written 20 months ago by _r_am32k

What do you mean by but these numbers mask the actual count numbers?

ADD REPLYlink written 20 months ago by ATpoint44k

Hi! I meant that after vsd transformation, genes with let say 1000 and 10 reads will be at the same scale. Nevertheless, the gene with 1000 reads was, probably, detected with higher accuracy, so I'd like to preserve such info. May be for this I should somehow switch to TPMs, but I'd like to avoid it for analysis transparence.


ADD REPLYlink written 20 months ago by Eugene A90

Hi Eugene,

how did you generate the box plot from the vsd data? thanks,

ADD REPLYlink written 11 months ago by afadda20

Looks like custom code using ggplot2

ADD REPLYlink written 11 months ago by ATpoint44k

something like this (given matrix with row - gnes, sample - columns):

 pritty_boxplot <- function(matrix, masker, row, deseq_res = FALSE){
          #print nice boxplot for given row in dataframe (matrix)
          #masker - columns lables

                 value <- as.numeric(as.vector(unlist(matrix[row,])))
                 value <- as.vector(unlist(counts(matrix[row], normalized = TRUE)))
         data <- data.frame(masker,value)
         data$lab <- colnames(matrix)
         ggplot(data, aes(y=value,x = masker ,fill=masker)) + geom_boxplot(outlier.shape = NA) +
                    geom_point(position = position_jitter(w = 0.2,seed = 2))+ 

         #to label points on the plot uncomment
         #geom_text(aes(label=lab),position =position_jitter(w = 0.2,seed = 2) ,angle=30) + 

         #scale_y_log10() + 
         xlab('') + ylab(row)
ADD REPLYlink modified 9 months ago • written 9 months ago by Eugene A90

Hi! I want to know how you made this boxplot, I saw the script that you used for that but a couldn't get some things. I guess that you converted vst object into matrix called "matrix", but what means "masker" and "row" and also that thing of "deseq_res = FALSE" in the script? I would appreciate if you could explain this for me, I am new in this! Thanks a lot! n_n

ADD REPLYlink written 12 weeks ago by fanimata0

matrix - is either count matrix or Deseq object - if later is true - set deseq_res = True

row - gene name you whant to plot 'RAG1'

masker - vector of the length == number of columns in your matrix, which contain group names you whant to see on a chart - c('contr', 'contr', 'contr', 'exp', 'exp') for example (in given example your matrix might have the following colnames: c('control_1', 'control_2', 'control_3', 'experiment_1', 'experiment_2')

ADD REPLYlink modified 12 weeks ago • written 12 weeks ago by Eugene A90
gravatar for ATpoint
20 months ago by
ATpoint44k wrote:

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.

ADD COMMENTlink modified 20 months ago • written 20 months ago by ATpoint44k

Thanks for the answer!! Shame but it did not occur to me that vst is indeed that much close to the log transformation. Best, Eugene

PS I'd mark your comment as an answer, if I could

ADD REPLYlink modified 20 months ago • written 20 months ago by Eugene A90
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1669 users visited in the last hour