DESeq2: Difference in rlog and vst values if run transform before or after running DESeq
Entering edit mode
Zach Roe ▴ 10
Last seen 18 months ago


I notice that the rlog and vst values for my assay differ depending if I run the DESeq function before or after running the transformation.

I thought that did this not matter as the manual states under section Data transformations and visualization -- "However, the running times are shorter and more similar with blind=FALSE and if the function DESeq has already been run, because then it is not necessary to re-estimate the dispersion values," so I assume the only difference was running time, and the transforms are performing the estimation of dispersion values if DESeq was not run prior?

Should DESeq be run before transforms?

# Case 1

dds1 <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~ condition)
rld1 <- rlog(dds1, blind=FALSE)
vsd1 <- varianceStabilizingTransformation(dds1, blind=FALSE)


# Case 2

dds2 <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~ condition)
dds2 <- DESeq(dds2)
rld2 <- rlog(dds2, blind=FALSE)
vsd2 <- varianceStabilizingTransformation(dds2, blind=FALSE)

length(which((assay(dds1) == assay(dd2)) == FALSE))
# [1] 0
length(which((assay(rld1) == assay(rld2)) == FALSE))
# [1] 583555 # the entire set
length(which((assay(vsd1) == assay(vsd2)) == FALSE))
# [1] 583555 # the entire set

I am running DESeq2_1.14.1 

R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.3

deseq2 rlog variancestabilizingtransformation DESeq • 1.7k views
Entering edit mode
Last seen 2 days ago
United States

Not sure why you're not getting equality, at least within tolerance. Can you use all.equal(x,y) instead, so we can see what the mean relative difference is?

Here's what I get:

> dds <- makeExampleDESeqDataSet()
> rld <- rlog(dds, blind=FALSE)
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> rld2 <- rlog(dds, blind=FALSE)
> all.equal(assay(rld), assay(rld2))

[1] TRUE

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.10

 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] DESeq2_1.14.1              SummarizedExperiment_1.4.0 Biobase_2.34.0            
 [4] GenomicRanges_1.26.1       GenomeInfoDb_1.10.1        IRanges_2.8.1             
 [7] S4Vectors_0.12.0           BiocGenerics_0.20.0        rmarkdown_1.2             
[10] magrittr_1.5               knitr_1.15.1               devtools_1.12.0           
[13] BiocInstaller_1.24.0      

loaded via a namespace (and not attached):
 [1] genefilter_1.56.0    locfit_1.5-9.1       splines_3.3.2        lattice_0.20-34     
 [5] colorspace_1.3-1     htmltools_0.3.5      survival_2.40-1      XML_3.98-1.5        
 [9] foreign_0.8-67       withr_1.0.2          DBI_0.5-1            BiocParallel_1.8.1  
[13] RColorBrewer_1.1-2   plyr_1.8.4           stringr_1.1.0        zlibbioc_1.20.0     
[17] munsell_0.4.3        gtable_0.2.0         memoise_1.0.0        evaluate_0.10       
[21] latticeExtra_0.6-28  geneplotter_1.52.0   AnnotationDbi_1.36.0 htmlTable_1.7       
[25] Rcpp_0.12.8          acepack_1.4.1        xtable_1.8-2         scales_0.4.1        
[29] backports_1.0.4      Hmisc_4.0-0          annotate_1.52.0      XVector_0.14.0      
[33] gridExtra_2.2.1      ggplot2_2.2.0        digest_0.6.10        stringi_1.1.2       
[37] grid_3.3.2           rprojroot_1.1        tools_3.3.2          bitops_1.0-6        
[41] lazyeval_0.2.0       RCurl_1.95-4.8       tibble_1.2           RSQLite_1.1         
[45] Formula_1.2-1        cluster_2.0.5        Matrix_1.2-8         data.table_1.9.8    
[49] assertthat_0.1       rpart_4.1-10         nnet_7.3-12 
Entering edit mode

Hi Michael I am still trying to figure out where the discrepancy occurs.  I had also been filtering using fpm and thought that was causing the error but I cannot reproduce the error when using the example DESeqDataSet.


Login before adding your answer.

Traffic: 312 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6