Limma: topTable without eBayes?
1
0
Entering edit mode
Guido Hooiveld ★ 3.9k
@guido-hooiveld-2020
Last seen 1 day ago
Wageningen University, Wageningen, the …
Dear list, topTable is a convenient function for displaying the top regulated genes in an experiment. topTable is usually applied on a object that is generated after fitting a linear model AND subsequent application of the eBayes funtion (fit3 below). I would like to view the results of a model WITHOUT applying the eBayes function, preferably using topTable. That is, using topTable to summerize the results of fit2 below. Applying topTable as such won't work; apparently a format imposed by the eBayes function is required to have topTable properly work. Even specifying a contrast of interest doesn't work. Therefore, has anyone a suggestion how to summerize the results of the fit2 mentioned below? Thanks, Guido > library(affy) > library(limma) > targets <- readTargets("targets_A23.txt") > > affy.data <- ReadAffy(filenames=targets$FileName) > x.norm <- rma(affy.data) Background correcting Normalizing Calculating Expression > > TS <- paste(targets$Strain, targets$Treatment, sep=".") > TS <- factor(TS, levels=c("WT.Con","WT.WY","KO.Con","KO.WY")) > design <- model.matrix(~0+TS) > colnames(design) <- levels(TS) > fit <- lmFit(x.norm, design) > cont.matrix <- makeContrasts(A=WT.WY-WT.Con, B=KO.WY-KO.Con, levels=design) > > fit2 <- contrasts.fit(fit, cont.matrix) #<<---- how to conviently summerize results of this fit? > > fit3 <- eBayes(fit2) > > topTable(fit3) ID A B AveExpr F P.Value adj.P.Val 13289 1431833_a_at 3.538370 -0.091190845 7.856618 1453.7786 5.084771e-18 1.153734e-13 18706 1450643_s_at 2.847376 0.182291682 7.907519 1234.5101 1.739465e-17 1.620566e-13 17177 1449065_at 4.528581 0.125074564 6.713129 1200.7506 2.142661e-17 1.620566e-13 7120 1422925_s_at 3.149781 -0.007867198 8.072523 895.1004 1.945926e-16 1.103827e-12 7192 1422997_s_at 4.298410 0.090012125 8.311977 837.9040 3.193016e-16 1.448991e-12 13047 1431012_a_at 2.682180 0.228653000 8.911287 785.9448 5.159295e-16 1.698832e-12 6721 1422526_at 2.772829 0.123117132 8.060643 784.2985 5.240997e-16 1.698832e-12 16494 1448382_at 3.386068 0.193350865 9.658716 685.0981 1.442685e-15 3.703026e-12 296 1415965_at 3.749143 0.071163645 5.293939 683.4572 1.468807e-15 3.703026e-12 8911 1424716_at 2.647513 0.370083995 8.453969 665.2223 1.798165e-15 4.080035e-12 > topTable(fit2) Error in topTableF(fit, number = number, genelist = genelist, adjust.method = adjust.method, : F-statistics not available. Try topTable for individual coef instead. > topTable(fit2, coef=1) Error in dim(data) <- dim : attempt to set an attribute on NULL > > sessionInfo() R version 2.11.0 (2010-04-22) i386-pc-mingw32 locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] moe430acdf_2.6.0 limma_3.4.0 affy_1.26.1 Biobase_2.8.0 loaded via a namespace (and not attached): [1] affyio_1.16.0 preprocessCore_1.10.0 tools_2.11.0 > ------------------------------------------------ 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<http: nutrigene.4t.com=""/> email: guido.hooiveld@wur.nl [[alternative HTML version deleted]]
affy affy • 2.1k views
ADD COMMENT
0
Entering edit mode
Jenny Drnevich ★ 2.0k
@jenny-drnevich-2812
Last seen 12 hours ago
United States
Hi Guido, topTable for a single coef needs the $t, $p.value and $lods that are created from eBayes. You could manually compute these to be regular t-statistics instead of those from eBayes. See section 10 of the limmaUsersGuide(). Try this (I didn't test it): fit2$t <- fit2$coef/fit2$stdev.unscaled/fit2$sigma fit2$p.value <- 2 * pt(-abs(fit2$t), df = fit2$df.residual) The $lods is harder; if you care about it you'll have to figure out how to duplicate it from the ebayes code by using modified inputs, which is what I did for the $t and $p-value; if you don't care about it but need something for topTable to work, I'd just do: fit2$lods <- fit3$lods HTH, Jenny At 12:51 PM 8/31/2010, Hooiveld, Guido wrote: >Dear list, > >topTable is a convenient function for displaying the top regulated >genes in an experiment. topTable is usually applied on a object that >is generated after fitting a linear model AND subsequent application >of the eBayes funtion (fit3 below). >I would like to view the results of a model WITHOUT applying the >eBayes function, preferably using topTable. That is, using topTable >to summerize the results of fit2 below. Applying topTable as such >won't work; apparently a format imposed by the eBayes function is >required to have topTable properly work. Even specifying a contrast >of interest doesn't work. >Therefore, has anyone a suggestion how to summerize the results of >the fit2 mentioned below? > >Thanks, >Guido > > > library(affy) > > library(limma) > > targets <- readTargets("targets_A23.txt") > > > > affy.data <- ReadAffy(filenames=targets$FileName) > > x.norm <- rma(affy.data) >Background correcting >Normalizing >Calculating Expression > > > > TS <- paste(targets$Strain, targets$Treatment, sep=".") > > TS <- factor(TS, levels=c("WT.Con","WT.WY","KO.Con","KO.WY")) > > design <- model.matrix(~0+TS) > > colnames(design) <- levels(TS) > > fit <- lmFit(x.norm, design) > > cont.matrix <- makeContrasts(A=WT.WY-WT.Con, B=KO.WY-KO.Con, levels=design) > > > > fit2 <- contrasts.fit(fit, cont.matrix) #<<---- how to > conviently summerize results of this fit? > > > > fit3 <- eBayes(fit2) > > > > topTable(fit3) > ID A B AveExpr F > P.Value adj.P.Val >13289 1431833_a_at 3.538370 -0.091190845 7.856618 1453.7786 >5.084771e-18 1.153734e-13 >18706 1450643_s_at 2.847376 0.182291682 7.907519 1234.5101 >1.739465e-17 1.620566e-13 >17177 1449065_at 4.528581 0.125074564 6.713129 1200.7506 >2.142661e-17 1.620566e-13 >7120 1422925_s_at 3.149781 -0.007867198 8.072523 895.1004 >1.945926e-16 1.103827e-12 >7192 1422997_s_at 4.298410 0.090012125 8.311977 837.9040 >3.193016e-16 1.448991e-12 >13047 1431012_a_at 2.682180 0.228653000 8.911287 785.9448 >5.159295e-16 1.698832e-12 >6721 1422526_at 2.772829 0.123117132 8.060643 784.2985 >5.240997e-16 1.698832e-12 >16494 1448382_at 3.386068 0.193350865 9.658716 685.0981 >1.442685e-15 3.703026e-12 >296 1415965_at 3.749143 0.071163645 5.293939 683.4572 >1.468807e-15 3.703026e-12 >8911 1424716_at 2.647513 0.370083995 8.453969 665.2223 >1.798165e-15 4.080035e-12 > > topTable(fit2) >Error in topTableF(fit, number = number, genelist = genelist, >adjust.method = adjust.method, : > F-statistics not available. Try topTable for individual coef instead. > > topTable(fit2, coef=1) >Error in dim(data) <- dim : attempt to set an attribute on NULL > > > > sessionInfo() >R version 2.11.0 (2010-04-22) >i386-pc-mingw32 > >locale: >[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United >States.1252 LC_MONETARY=English_United States.1252 >LC_NUMERIC=C LC_TIME=English_United States.1252 > >attached base packages: >[1] stats graphics grDevices utils datasets methods base > >other attached packages: >[1] moe430acdf_2.6.0 limma_3.4.0 affy_1.26.1 Biobase_2.8.0 > >loaded via a namespace (and not attached): >[1] affyio_1.16.0 preprocessCore_1.10.0 tools_2.11.0 > > > > >------------------------------------------------ >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<http: nutrigene.4t.com=""/> >email: guido.hooiveld at wur.nl > > > > [[alternative HTML version deleted]] > >_______________________________________________ >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
ADD COMMENT

Login before adding your answer.

Traffic: 554 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6