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