edgeR - one way ANOVA, coef parameter interpretation in glmLRT()
1
0
Entering edit mode
@michal-okoniewski-2676
Last seen 10.2 years ago
Hey Mark, A question on one-way ANOVA design with edgeR: I do something like that for 3 groups, so similar case like in the vignette (glm functionality): d <- calcNormFactors(d) d <- estimateCommonDisp(d, verbose=TRUE) d <- estimateTagwiseDisp(d) TS <- factor(cc) design<- model.matrix(~0+TS) colnames(design) <- levels(TS) fit <- glmFit(d, design) fit2 <- glmLRT(fit) # or fit2 <- glmLRT(fit, coef=1) ??? outTab <- topTags(fit, n=30666) An the question is: __ How to get the p-values of the F-test for one-way anova? In a basic comparison of 3 groups? __ In the vignette you have written: "The fi t has three parameters. The fi rst is the baseline level of group 1. The second and third are the 2 vs 1 and 3 vs 1 di ferences." edgeR vignette (p17, date March31st 2013) This sounds to me a bit like hard-coded contrasts? in addition, default value of coef is ncol(glmfit$design), so 3 in this case. Shall I use coef=1? I am a bit confused? Thanks! I send a copy to BioC list, as I could not find the answer there. If I ask stupid question ? please give me a simple solution too :) Cheers, Michal > sessionInfo() R version 3.0.1 (2013-05-16) Platform: x86_64-apple-darwin10.8.0 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] edgeR_3.2.3 limma_3.16.5 loaded via a namespace (and not attached): [1] tools_3.0.1 -- Dr. Michal Okoniewski Functional Genomics Center Zurich (mostly Wed-Fri) Winterthurerstrasse 190, 8057 Zurich Room Y32 H66, Phone: +41 44 635 39 24 Department of Neuroimmunology and MS Research, UniSpital Zurich, (mostly Mon-Thu) Raemistrasse 100, 8091 Zurich Room: OPS D37, Phone: +41 44 255 1756
edgeR edgeR • 2.1k views
ADD COMMENT
0
Entering edit mode
Mark Robinson ▴ 880
@mark-robinson-4908
Last seen 6.1 years ago
Hey Michal, I've added some comments in line below ? On 10.06.2013, at 15:34, Michal Okoniewski <michal.okoniewski at="" fgcz.ethz.ch=""> wrote: > Hey Mark, > > A question on one-way ANOVA design with edgeR: > > I do something like that for 3 groups, so similar case like in the vignette (glm functionality): > > > d <- calcNormFactors(d) > > d <- estimateCommonDisp(d, verbose=TRUE) > > d <- estimateTagwiseDisp(d) > > > TS <- factor(cc) > > design<- model.matrix(~0+TS) > > colnames(design) <- levels(TS) You seem to have classic edgeR (from 2.7, [1]) but probably want GLM edgeR (from 2.8, [1]), no? See below. > fit <- glmFit(d, design) > > fit2 <- glmLRT(fit) > > # or fit2 <- glmLRT(fit, coef=1) ??? > > > outTab <- topTags(fit, n=30666) > > An the question is: > __ How to get the p-values of the F-test for one-way anova? In a basic comparison of 3 groups? __ > > In the vignette you have written: > > "The fi t has three parameters. The fi rst is the baseline level of group 1. The second and third > are the 2 vs 1 and 3 vs 1 di ferences." edgeR vignette (p17, date March31st 2013) > > This sounds to me a bit like hard-coded contrasts? in addition, default value of coef is ncol(glmfit$design), so 3 in this case. Shall I use coef=1? > I am a bit confused? > > Thanks! You need to be careful of what you use for the design matrix, or more precisely, how you pair your design matrix with your contrast matrix (or what coefficients of the design matrix you specify). Let me try and spell it out. In the vignette (2.8.3, [1]), the design matrix is defined as: group <- factor(c(1,1,2,2,3,3)) design <- model.matrix(~group) > design (Intercept) group2 group3 1 1 0 0 2 1 0 0 3 1 1 0 4 1 1 0 5 1 0 1 6 1 0 1 Note that your design matrix, as defined above, is *different*, like: group <- factor(c(1,1,2,2,3,3)) design <- model.matrix(~0+group) > design group1 group2 group3 1 1 0 0 2 1 0 0 3 0 1 0 4 0 1 0 5 0 0 1 6 0 0 1 So, what to do next (in terms of glmLRT) depends on how you've made the design matrix. For a 1-way ANOVA style comparison with 3 levels of the factor, I find the following the simplest: d <- DGEList(...) d <- calcNormFactors(d) ----- design <- model.matrix(~group) d <- estimateGLMCommonDisp(d, design) d <- estimateGLMTrendedDisp(d, design) d <- estimateGLMTagwiseDisp(d, design) f <- glmFit(d, design) lrt <- glmLRT(f, coef=2:3) topTags(lrt) ----- but note that this is equivalent to: ----- design <- model.matrix(~0+group) d <- estimateGLMCommonDisp(d, design) d <- estimateGLMTrendedDisp(d, design) d <- estimateGLMTagwiseDisp(d, design) f <- glmFit(d, design) lrt <- glmLRT(f, contrast=matrix(c(-1,1,0,0,-1,1),ncol=2)) topTags(lrt) ----- So, really, it's a matter of pairing your design matrix with the correct 'coef' or 'contrast' argument of glmLRT(). Hope that helps, Mark [1] http://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/ inst/doc/edgeRUsersGuide.pdf > I send a copy to BioC list, as I could not find the answer there. > If I ask stupid question ? please give me a simple solution too :) > > Cheers, > Michal > > >> sessionInfo() > R version 3.0.1 (2013-05-16) > Platform: x86_64-apple-darwin10.8.0 (64-bit) > > locale: > [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] edgeR_3.2.3 limma_3.16.5 > > loaded via a namespace (and not attached): > [1] tools_3.0.1 > > > > -- > Dr. Michal Okoniewski > > Functional Genomics Center Zurich (mostly Wed-Fri) > Winterthurerstrasse 190, 8057 Zurich > Room Y32 H66, Phone: +41 44 635 39 24 > > Department of Neuroimmunology and MS Research, > UniSpital Zurich, (mostly Mon-Thu) > Raemistrasse 100, 8091 Zurich > Room: OPS D37, Phone: +41 44 255 1756
ADD COMMENT

Login before adding your answer.

Traffic: 724 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