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

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).

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.