[BioC] a possible bug in DESeq

[Wolfgang Huber]

Shan Gao

you have sent me your file "expression-qwei.txt". The code that you 
posted in your initial mail could not possibly have worked with it, 
since it contains 3 columns, while your 'conditions' vector only has 
length 2. I have run the below script (which approximates as much as I 
could tell what you claimed you did) without problem. Its output is also 
posted below, as is the session info.

I conclude that this was operator error, and there is no such bug in DESeq.

I noted that your session info mentions a version of DESeq (1.6.1) that 
does not match your version of R (2.15.1). It seems that your R 
installation is messed up. Perhaps you can ask your system administrator 
to give you a correct installation of R and Bioconductor.


Code:
-----
library("DESeq")

raw.data <- read.table("expression-qwei.txt", row.names=1)

conditions = c("mutant","wild")
counts = raw.data[, conditions]
counts = counts[rowSums(counts) >= length(conditions)/2,]
cds  = newCountDataSet(counts, conditions)

cds = estimateSizeFactors(cds)
print(sizeFactors(cds))

cds = estimateDispersions(cds, method='blind')
print(dispTable(cds))

cds = estimateDispersions(cds, method='blind', sharingMode='fit-only')
print(dispTable(cds))

Output:
-------

sessionInfo()
   mutant     wild
1.075829 0.929516
  mutant    wild
"blind" "blind"
  mutant    wild
"blind" "blind"

 > sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-apple-darwin11.4.0/x86_64 (64-bit)

locale:
[1] C

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

other attached packages:
[1] DESeq_1.8.3        locfit_1.5-8       Biobase_2.16.0 
BiocGenerics_0.2.0
[5] fortunes_1.5-0

loaded via a namespace (and not attached):
  [1] AnnotationDbi_1.18.1 DBI_0.2-5            IRanges_1.14.4
  [4] RColorBrewer_1.0-5   RSQLite_0.11.1       XML_3.9-4
  [7] annotate_1.34.1      genefilter_1.38.0    geneplotter_1.34.0
[10] grid_2.15.1          lattice_0.20-6       splines_2.15.1
[13] stats4_2.15.1        survival_2.36-14     xtable_1.7-0

On 7/29/12 7:30 PM, Wolfgang Huber wrote:
> Dear Shan Gao
>
> thank you. I was not able to reproduce what you report, based on the
> code below and a different data set. Can you please send us the data
> file "expression-qwei.txt" (or any other data set that for you produces
> this problem) together with the exact code that you ran (e.g. according
> to the output of the R function 'history') ?
>
>      Best wishes
>      Wolfgang
>
>
> On 7/29/12 5:13 PM, wang peter wrote:
>> i have two sample, one is wild, the other is mutant
>> i met a very wired problem
>> when i set parameter sharingMode = "fit-only", such coding cannot work
>> but if i set sharingMode=default, it can work
>>
>>
>> rm(list=ls())
>> library(DESeq)
>>
>> raw.data <- read.table("expression-qwei.txt",row.names=1) #
>> counts <- raw.data[, 1:dim(raw.data)[2]]#
>> conditions=c("mutant","wild") #
>> counts <- counts[rowSums(counts) >= length(conditions)/2,]#
>> cds <- newCountDataSet(counts, conditions)#
>>
>> #normalization
>> cds <- estimateSizeFactors(cds)
>> sizeFactors(cds)
>> cds <- estimateDispersions(cds, method = "blind"，sharingMode =
>> "fit-only")
>>> dispTable(cds)
>> [1] NA NA
>>
>>> cds <- estimateDispersions(cds, method = "blind")
>>> dispTable(cds)
>>   mutant    wild
>> "blind" "blind"
>>
>>> sessionInfo()
>> R version 2.15.1 (2012-06-22)
>> Platform: x86_64-pc-mingw32/x64 (64-bit)
>>
>> locale:
>> [1] LC_COLLATE=Chinese (Simplified)_People's Republic of China.936
>> LC_CTYPE=Chinese (Simplified)_People's Republic of China.936
>> [3] LC_MONETARY=Chinese (Simplified)_People's Republic of China.936
>> LC_NUMERIC=C
>> [5] LC_TIME=Chinese (Simplified)_People's Republic of China.936
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> other attached packages:
>> [1] DESeq_1.6.1    locfit_1.5-8   Biobase_2.14.0
>>
>> loaded via a namespace (and not attached):
>>   [1] annotate_1.32.3       AnnotationDbi_1.16.19 DBI_0.2-5
>> genefilter_1.36.0     geneplotter_1.32.1    grid_2.15.1
>> IRanges_1.12.6
>>   [8] lattice_0.20-6        RColorBrewer_1.0-5    RSQLite_0.11.1
>> splines_2.15.1        survival_2.36-14      xtable_1.7-0
>>