Hi all

Idiocy alert, I think...!

I am running a very simple differential expression experiment.  There are two conditions, neg and pos with three replicates of each and the experimental design is a simple ~ condition.  I am using counts from htseq-count.

For whatever reason, DESeq2 is not using a beta prior and the log2 fold changes generated are unshrunken - I am not sure why this is happening.

This is my code:

library("DESeq2")

sampleFiles <- list.files(pattern="\\.txt$")

sampleCondition <- c("neg", "neg", "pos", "pos",  "neg", "pos")
sampleTable <- data.frame(sampleName = sampleFiles, fileName = sampleFiles, condition=sampleCondition)

ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory=getwd(), design = ~ condition)
ddsHTSeq <- ddsHTSeq[ rowSums(counts(ddsHTSeq)) >1, ]

ddsHTSeq$condition <- relevel(ddsHTSeq$condition, ref="neg")

dds <- DESeq(ddsHTSeq)
res <- results(dds, alpha=0.5)

plotMA(res)

The sampleTable looks like this:

> sampleTable
                         sampleName                          fileName condition
1 17970-24hneg2-34062058_counts.txt 17970-24hneg2-34062058_counts.txt       neg
2 17970-24hneg3-34037178_counts.txt 17970-24hneg3-34037178_counts.txt       neg
3 17970-24hpos2-34036147_counts.txt 17970-24hpos2-34036147_counts.txt       pos
4 17970-24hpos3-34062062_counts.txt 17970-24hpos3-34062062_counts.txt       pos
5 19047-24hneg1-34050175_counts.txt 19047-24hneg1-34050175_counts.txt       neg
6 19047-24hpos1-34039162_counts.txt 19047-24hpos1-34039162_counts.txt       pos

The MA plot looks like this: http://imgur.com/a/1rNLh (sorry, couldn't get this to embed!)

If I try to addMLE I get an error implying that the beta prior is not being used:

res <- results(dds, alpha=0.05, addMLE = TRUE)

Error in results(dds, alpha = 0.05, addMLE = TRUE) : 
  addMLE=TRUE is only for when a beta prior was used.
otherwise, the log2 fold changes are already MLE

I have analysed this dataset before with no problems, but now I am doing a re-analysis to incorporate the third replicate and I am not sure why this is happening.  I must be doing something idiotic - any insight appreciated!

Thanks!

> sessionInfo()
R version 3.4.1 (2017-06-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.5 LTS

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.0
LAPACK: /usr/lib/lapack/liblapack.so.3.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] DESeq2_1.16.1              SummarizedExperiment_1.6.3 DelayedArray_0.2.7         matrixStats_0.52.2         Biobase_2.36.2             GenomicRanges_1.28.4      
 [7] GenomeInfoDb_1.12.2        IRanges_2.10.2             S4Vectors_0.14.3           BiocGenerics_0.22.0       

loaded via a namespace (and not attached):
 [1] genefilter_1.58.1       locfit_1.5-9.1          splines_3.4.1           lattice_0.20-35         colorspace_1.3-2        htmltools_0.3.6        
 [7] base64enc_0.1-3         blob_1.1.0              survival_2.41-3         XML_3.98-1.9            rlang_0.1.2             DBI_0.7                
[13] foreign_0.8-69          BiocParallel_1.10.1     bit64_0.9-7             RColorBrewer_1.1-2      GenomeInfoDbData_0.99.0 plyr_1.8.4             
[19] stringr_1.2.0           zlibbioc_1.22.0         munsell_0.4.3           gtable_0.2.0            htmlwidgets_0.9         memoise_1.1.0          
[25] latticeExtra_0.6-28     knitr_1.17              geneplotter_1.54.0      AnnotationDbi_1.38.2    htmlTable_1.9           Rcpp_0.12.12           
[31] acepack_1.4.1           xtable_1.8-2            scales_0.4.1            backports_1.1.0         checkmate_1.8.3         Hmisc_4.0-3            
[37] annotate_1.54.0         XVector_0.16.0          bit_1.1-12              gridExtra_2.2.1         ggplot2_2.2.1           digest_0.6.12          
[43] stringi_1.1.5           grid_3.4.1              tools_3.4.1             bitops_1.0-6            magrittr_1.5            RSQLite_2.0            
[49] lazyeval_0.2.0          RCurl_1.95-4.8          tibble_1.3.3            Formula_1.2-2           cluster_2.0.6           Matrix_1.2-10          
[55] data.table_1.10.4       rpart_4.1-11            nnet_7.3-12             compiler_3.4.1         

Take a look at the NEWS file for DESeq2. The shrinkage step has been pulled out of the analysis pipeline, and can be applied to your result to perform shrinkage.

Read through the vignette again, particularly hunting for the use of the lfcShrink funciton.

I believe Mike has previously mentioned elsewhere on this forum that at least part of his motivation for this was to set the stage to make it easier to experiment with different types of shrinkage ...