[BioC] Testing main effects and interactions in edgeR without contrasts
David Shuker
dms14 at st-andrews.ac.uk
Sat Dec 14 18:57:58 CET 2013
Dear Gordon,
Many thanks for your email, and also to Jim for helping to clarify the syntax for removing terms.
Very much appreciated.
Best,
Dave
Dear David,
Here's a toy example showing how to use edgeR to test for interactions and
main effects. I'll generate Poisson counts so we can compare with glm()
in the stats package. Generate counts for two genes and two replicates
per treatment combination:
y <- matrix(rpois(24,lambda=20),2,12)
A <- gl(2,3,12)
B <- gl(3,1,12)
Assume library sizes are all equal to 10,000.
lib.size <- rep(10000,12)
offset <- log(lib.size)
The standard factorial analysis for the first gene would be as follows:
out <- glm(y[1,] ~ A*B, family="poisson", offset=offset)
anova(out, test="Chi")
When I ran the above code, it gave me the following Analysis of Deviance
Table:
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 11 25.538
A 1 0.7016 10 24.837 0.4023
B 2 1.4055 8 23.431 0.4952
A:B 2 8.6404 6 14.790 0.0133 *
For the second gene, I got:
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 11 9.2948
A 1 2.1341 10 7.1607 0.1441
B 2 1.4481 8 5.7127 0.4848
A:B 2 1.4377 6 4.2750 0.4873
Now test for interactions using edgeR:
design <- model.matrix(~A*B)
fit <- glmFit(y,design,offset=offset,dispersion=0)
lrt <- glmLRT(fit,coef=5:6)
topTags(lrt)
This gives:
Coefficient: A2:B2 A2:B3
logFC.A2.B2 logFC.A2.B3 logCPM LR PValue FDR
1 -0.78611 0.59857 11.108 8.6404 0.013297 0.026594
2 0.36661 -0.22905 10.910 1.4377 0.487313 0.487313
Notice that edgeR gives the identical chisquare statistics (8.6404 and
1.4377) and the identical p-values (0.0133 and 0.4873) as does anova() for
the interaction A:B.
To test for main effects, you would fit the additive model:
design <- model.matrix(~A+B)
fit <- glmFit(y,design,offset=offset,dispersion=0)
Then
topTags(glmLRT(fit,coef=3:4))
Coefficient: B2 B3
logFC.B2 logFC.B3 logCPM LR PValue FDR
1 -0.20769 -0.245411 11.108 1.4055 0.49521 0.49521
2 -0.23739 0.039259 10.910 1.4481 0.48479 0.49521
gives the exactly the same statistics for B as in the ANODEV tables.
Best wishes
Gordon
> Date: Fri, 13 Dec 2013 02:25:36 -0800 (PST)
> From: "David Shuker [guest]" <guest at bioconductor.org>
> To: bioconductor at r-project.org, david.shuker at st-andrews.ac.uk
> Subject: [BioC] Testing main effects and interactions in edgeR without
> contrasts
>
> Dear friends,
>
> We are very interested in testing main effects and interactions for a
> 2x3 factorial RNA-seq project we have been running (we have reasonable
> replication: 7 biological replicates per treatment combination, i.e. 42
> libraries in total).
>
>> From our reading of the edgeR Users guide and various posts on-line, it
>> looks as though the package is set-up to test effects via contrasts.
>> Although a contrasts-approach allows the testing of differences between
>> specific treatment combinations, we would like to test interactions and
>> main effects in a more traditional glm/anova format. (We note that
>> testing whether the two contrasts associated with a three-level main
>> effect - a, b, c - are significant is not the same question as testing
>> whether the main effect itself is significant, or of course vice
>> versa.)
>
> We would like outputs and tests along the lines of glm() or anova(), not
> "anova-like" tests that use the results of one or more contrasts to ask
> similar, but not exactly the same, questions of the data.
>
> We apologise in advance if we have either missed something obvious in
> the user guide or relevant discussion somewhere online. If edgeR does
> not have this functionality, we will probably use edgeR in conjunction
> with other packages.
>
> As with all the community, we remain in awe of the work the edgeR
> developers have done, and thank them in advance.
>
> Best wishes,
>
> Dave Shuker
>
> -- output of sessionInfo():
>
> NA
>
> --
> Sent via the guest posting facility at bioconductor.org.
______________________________________________________________________
The information in this email is confidential and intend...{{dropped:6}}
More information about the Bioconductor
mailing list