[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