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

Laurent Gautier laurent at cbs.dtu.dk
Sun Jun 14 10:03:59 CEST 2009


Hooiveld, Guido wrote:
> 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.

It does not seem so.

There is a default value (zero), both according to the man page and to 
signature of the function:

 > args(treat)
function (fit, lfc = 0)
NULL


> - you have to run the eBayes function on 'fit'.


Thanks.

If this is the case, the man page should probably specify it.
Currently there is:

"
treat(fit, lfc=0)

Arguments

fit	 an MArrayLM fitted model object produced by lmFit or contrasts.fit, 
or an unclassed list produced by lm.series, gls.series or mrlm 
containing components coefficients, stdev.unscaled, sigma and df.residual
"

Also, the man page says further down:

"
treat is concerned with p-values rather than posterior odds, so it does 
not compute the B-statistic lods.
"
Shouldn't the column B be discarded by treat() then ?




L.



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



More information about the Bioconductor mailing list