[DESeq2] Different results from tutorial
2
0
Entering edit mode
@changhan1110-15895
Last seen 5.9 years ago

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

deseq2 tximportdata tximport • 1.6k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode
ADD REPLY
2
Entering edit mode
@mikelove
Last seen 5 hours ago
United States

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 COMMENT
0
Entering edit mode

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 REPLY
0
Entering edit mode
@wolfgang-huber-3550
Last seen 3 months ago
EMBL European Molecular Biology Laborat…

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 COMMENT
0
Entering edit mode

Thanks for your comment. I will check it.

ADD REPLY

Login before adding your answer.

Traffic: 620 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

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

Powered by the version 2.3.6