Question: DESeq2: Difference in rlog and vst values if run transform before or after running DESeq
gravatar for Zach Roe
21 months ago by
Zach Roe10
Zach Roe10 wrote:


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

ADD COMMENTlink modified 21 months ago by Michael Love20k • written 21 months ago by Zach Roe10
gravatar for Michael Love
21 months ago by
Michael Love20k
United States
Michael Love20k wrote:

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 
ADD COMMENTlink written 21 months ago by Michael Love20k

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.

ADD REPLYlink modified 20 months ago • written 20 months ago by Zach Roe10
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.2.0
Traffic: 205 users visited in the last hour