Search
Question: [DESeq2] Different results from tutorial
0
gravatar for changhan1110
11 weeks ago by
changhan11100 wrote:

Hi all,

I'm stuck in a strange result while doing DESeq2 tutorial with tximport and tximportData. 

The script I used is below:

library("DESeq2")
library("tximport")
library("tximportData")
dir <- system.file("extdata", package = "tximportData")
list.files(dir)
samples <- read.table(file.path(dir, "samples.txt"), header = TRUE)
files <- file.path(dir, "rsem", samples$run, paste0(samples$run, ".genes.results.gz"))
names(files) <- paste0("sample", 1:6)
txi.rsem <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE)
txi.rsem$length[txi.rsem$length == 0] <- 1
sampleTable <- data.frame(condition = factor(rep(c("A", "B"), each = 3)))
rownames(sampleTable) <- colnames(txi.rsem$counts)
dds <- DESeqDataSetFromTximport(txi.rsem, sampleTable, ~ condition)
dds <- DESeq(dds)
res <- results(dds)
resOrdered <- res[order(res$padj),]
summary(res)

And the result is:

## out of 29939 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 3, 0.01%
## LFC < 0 (down)     : 0, 0%
## outliers [1]       : 82, 0.27%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

What I want is :

## 
## out of 9921 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 518, 5.2%
## LFC < 0 (down)     : 536, 5.4%
## outliers [1]       : 1, 0.01%
## low counts [2]     : 1539, 16%
## (mean count < 6)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

> sessionInfo()

R version 3.5.0 (2018-04-23)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=Korean_Korea.949  LC_CTYPE=Korean_Korea.949   
[3] LC_MONETARY=Korean_Korea.949 LC_NUMERIC=C                
[5] LC_TIME=Korean_Korea.949    

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

other attached packages:
 [1] tximportData_1.8.0          tximport_1.8.0             
 [3] DESeq2_1.20.0               stringi_1.1.7              
 [5] survival_2.41-3             SummarizedExperiment_1.10.1
 [7] DelayedArray_0.6.0          BiocParallel_1.14.1        
 [9] matrixStats_0.53.1          Biobase_2.40.0             
[11] GenomicRanges_1.32.3        GenomeInfoDb_1.16.0        
[13] IRanges_2.14.10             S4Vectors_0.18.2           
[15] BiocGenerics_0.26.0        

loaded via a namespace (and not attached):
 [1] bit64_0.9-7            splines_3.5.0         
 [3] Formula_1.2-3          latticeExtra_0.6-28   
 [5] blob_1.1.1             GenomeInfoDbData_1.1.0
 [7] yaml_2.1.19            pillar_1.2.2          
 [9] RSQLite_2.1.1          backports_1.1.2       
[11] lattice_0.20-35        digest_0.6.15         
[13] RColorBrewer_1.1-2     XVector_0.20.0        
[15] checkmate_1.8.5        colorspace_1.3-2      
[17] htmltools_0.3.6        Matrix_1.2-14         
[19] plyr_1.8.4             pkgconfig_2.0.1       
[21] XML_3.98-1.11          genefilter_1.62.0     
[23] zlibbioc_1.26.0        xtable_1.8-2          
[25] scales_0.5.0           htmlTable_1.11.2      
[27] tibble_1.4.2           annotate_1.58.0       
[29] ggplot2_2.2.1          nnet_7.3-12           
[31] lazyeval_0.2.1         magrittr_1.5          
[33] memoise_1.1.0          foreign_0.8-70        
[35] tools_3.5.0            data.table_1.11.2     
[37] hms_0.4.2              stringr_1.3.1         
[39] munsell_0.4.3          locfit_1.5-9.1        
[41] cluster_2.0.7-1        AnnotationDbi_1.42.1  
[43] compiler_3.5.0         rlang_0.2.0           
[45] grid_3.5.0             RCurl_1.95-4.10       
[47] rstudioapi_0.7         htmlwidgets_1.2       
[49] bitops_1.0-6           base64enc_0.1-3       
[51] gtable_0.2.0           DBI_1.0.0             
[53] R6_2.2.2               gridExtra_2.3         
[55] knitr_1.20             bit_1.1-13            
[57] Hmisc_4.1-1            readr_1.1.1           
[59] Rcpp_0.12.17           geneplotter_1.58.0    
[61] rpart_4.1-13           acepack_1.4.1

 

Are there any problems in my script?

Thanks.

Changhan

ADD COMMENTlink modified 11 weeks ago by Michael Love19k • written 11 weeks ago by changhan11100

When you say "What I want is", what is the basis for this desire?

ADD REPLYlink written 11 weeks ago by Wolfgang Huber13k

It is a result of an example from DESeq2 vignette. (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#p-values-and-adjusted-p-values)

ADD REPLYlink written 11 weeks ago by changhan11100
2
gravatar for Michael Love
11 weeks ago by
Michael Love19k
United States
Michael Love19k wrote:

These are different datasets (and organisms). The DESeq2 vignette uses pasilla (fly) while the tximport vignette shows how to import 6 samples from GEUVADIS (human).

ADD COMMENTlink written 11 weeks ago by Michael Love19k

Thanks for reply! I missed it. So, my result seems to have no problem at all with GEUVADIS sample? I got the same result using edgeR.

ADD REPLYlink written 11 weeks ago by changhan11100
0
gravatar for Wolfgang Huber
11 weeks ago by
EMBL European Molecular Biology Laboratory
Wolfgang Huber13k wrote:

Did you do already do what the output suggests? i.e.

[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

Since this seems to be one, if not the, cause of the difference between your outputs.

 

 

ADD COMMENTlink written 11 weeks ago by Wolfgang Huber13k

Thanks for your comment. I will check it.

ADD REPLYlink written 11 weeks ago by changhan11100
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: 145 users visited in the last hour