[BioC] Using limma for omnibus F tests

[Gavin Kelly]

Dear all,
I'm trying to find genes that fail the null hypothesis that they are all
the same across >2 treatment groups using the limma package (which I
find exceptionally useful, by the way).

If I were doing this for a single gene, I'd use a traditional F test;
limma's eBayes generates an F value (via classifyTestsF) from the
individual t statistics, but if I use a model without an intercept, I
believe the F reported will favour genes that are zero in all treatment
groups.  And if I use an intercept I'm struggling to make any contrasts
or manipulations of the design matrix that respect the required symmetry
between all treatment groups

I'd like a way to more generally report an F that reflects genes that
are the same (not necessarily zero) in all groups;  I suspect there's a
clever way to do this using contrasts.fit (I thought maybe looking at
all pairwise comparisons which doesn't work in the snippet belo) but
finding that is predicated on my being clever!  Anyone got any
suggestions?

library(limma)
set.seed(1)
pdat <- expand.grid(treatment=c("a","b","c"), rep=c("r1","r2","r3"))
sig <- matrix(0,nrow=50, ncol=nrow(pdat))
sig[1:25, pdat$treatment=="a"] <-  3 # 25 are up in a
sig[26:50, pdat$treatment=="b"] <- 3 # next 25 are up in b
err <- matrix(rnorm(25*nrow(pdat)), nrow=25, ncol=nrow(pdat))
dat <- sig + rbind(err, err)  #add the same noise to both 
dat <- rbind(dat, err, err, err, err) # add some non differentials

mdl <- model.matrix(~treatment-1, pdat)
ctr <- makeContrasts(treatmentb-treatmenta, treatmentc-treatmenta,
treatmentc-treatmentb,levels=mdl)
fit <- eBayes(contrasts.fit(lmFit(dat, mdl), ctr))
plot(fit$F) # I'd like this to repeat 1:25 == 26:50

Regards - Gavin
-- 
Gavin Kelly
Bioinformatics & Biostatistics
Cancer Research UK