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:4}}
