[BioC] Error when calling LIMMA's topTable() with an object returned by treat()

[Hooiveld, Guido]

Hi Laurent,
I recently used limma's TREAT argument without any problems. 
After comparing your code with mine I do have 2 remarks:
- you have to specify a cut-off value for treat.
- you have to run the eBayes function on 'fit'.

HTH,
Guido


Taking your example:

> library(limma)
> sd <- 0.3*sqrt(4/rchisq(100,df=4))
> y <- matrix(rnorm(100*6,sd=sd),100,6)
> rownames(y) <- paste("Gene",1:100)
> y[1:2,4:6] <- y[1:2,4:6] + 2
> design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
> options(digit=3)
> 
> fit <- lmFit(y,design)
>
> #apply the eBayes function on fit 
> fit2 <- eBayes(fit)
>
> #apply the treat function with argument a FC cut-off of 1.1 
> trfit <- treat(fit2, lfc=log2(1.1))
> 
> topTable(trfit, coef=2, adjust="BH")
          ID      logFC         t      P.Value  adj.P.Val           B
2     Gene 2  1.9577409  5.227446 0.0007778008 0.07778008 -0.02011076
1     Gene 1  1.2639445  4.147762 0.0036094095 0.18047048 -1.44624165
100 Gene 100  0.8883699  3.451312 0.0111156128 0.37052043 -2.45955109
62   Gene 62  0.6827269  3.344943 0.0154085696 0.38521424 -2.61953235
69   Gene 69  0.7037365  2.732057 0.0342969691 0.68593938 -3.55708668
11   Gene 11 -0.5883241 -2.431874 0.0570737751 0.71116024 -4.01793218
33   Gene 33 -0.5164955 -2.275967 0.0757425135 0.71116024 -4.25471971
98   Gene 98 -0.8803174 -2.271460 0.0596925864 0.71116024 -4.26152123
60   Gene 60  0.4559045  2.253581 0.0853392282 0.71116024 -4.28847799
66   Gene 66  0.5776649  2.251627 0.0728599645 0.71116024 -4.29142209
> topTable(fit2, coef=2, adjust="BH")
          ID      logFC         t      P.Value  adj.P.Val           B
2     Gene 2  1.9577409  5.227446 0.0006950502 0.06950502 -0.02011076
1     Gene 1  1.2639445  4.147762 0.0029382872 0.14691436 -1.44624165
100 Gene 100  0.8883699  3.451312 0.0081400748 0.23911671 -2.45955109
62   Gene 62  0.6827269  3.344943 0.0095646683 0.23911671 -2.61953235
69   Gene 69  0.7037365  2.732057 0.0247864409 0.47441398 -3.55708668
11   Gene 11 -0.5883241 -2.431874 0.0398870275 0.47441398 -4.01793218
33   Gene 33 -0.5164955 -2.275967 0.0510964049 0.47441398 -4.25471971
98   Gene 98 -0.8803174 -2.271460 0.0514632476 0.47441398 -4.26152123
60   Gene 60  0.4559045  2.253581 0.0529445231 0.47441398 -4.28847799
66   Gene 66  0.5776649  2.251627 0.0531089861 0.47441398 -4.29142209
> 
By comparing the p-values between fit2 and trfit you clearly observe the
effect of TREAT; with TREAT the p-values slightly higher (=less
significant).

 
------------------------------------------------ 
Guido Hooiveld, PhD 
Nutrition, Metabolism & Genomics Group
Division of Human Nutrition 
Wageningen University 
Biotechnion, Bomenweg 2 
NL-6703 HD Wageningen 
the Netherlands 
tel: (+)31 317 485788 
fax: (+)31 317 483342 
internet:   http://nutrigene.4t.com
email:      guido.hooiveld at wur.nl



> -----Original Message-----
> From: bioconductor-bounces at stat.math.ethz.ch 
> [mailto:bioconductor-bounces at stat.math.ethz.ch] On Behalf Of 
> Laurent Gautier
> Sent: 13 June 2009 13:05
> To: bioconductor mail list bioconductor at stat.math.ethz.ch
> Subject: [BioC] Error when calling LIMMA's topTable() with an 
> object returned by treat()
> 
> Dear list,
> 
> 
> Calling LIMMA's topTable() function with an object returned by
> treat() generates an error:
> 
> Error in dim(data) <- dim : attempt to set an attribute on NULL
> 
> 
> # example
> 
> sd <- 0.3*sqrt(4/rchisq(100,df=4))
> y <- matrix(rnorm(100*6,sd=sd),100,6)
> rownames(y) <- paste("Gene",1:100)
> y[1:2,4:6] <- y[1:2,4:6] + 2
> design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))
> options(digit=3)
> 
> fit <- lmFit(y,design)
> trfit <- treat(fit)
> topTable(trfit,coef=2)
> 
> 
> Does anyone have a workaround ?
> 
> 
> 
> Laurent
> 
> 
>  > sessionInfo()
> R version 2.9.0 (2009-04-17)
> i386-apple-darwin8.11.1
> 
> locale:
> C
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> 
> other attached packages:
> [1] Biostrings_2.12.1    IRanges_1.2.1        lattice_0.17-22 
> preprocessCore_1.6.0 limma_2.18.0         Biobase_2.4.1
> 
> loaded via a namespace (and not attached):
> [1] grid_2.9.0
>  >
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: 
> http://news.gmane.org/gmane.science.biology.informatics.conductor
> 
>