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

Hi,

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 8 months ago by Michael Love14k • written 8 months ago by Zach Roe10
1
gravatar for Michael Love
8 months ago by
Michael Love14k
United States
Michael Love14k 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

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=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 8 months ago by Michael Love14k

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 8 months ago • written 8 months ago by Zach Roe10
Please log in to add an answer.

Help
Access

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