9 days ago by

Cambridge, United Kingdom

Use `topTable`

to do your F-tests, don't try to poke around in the internals of the `MArrayLM`

object unless you know what you're doing. The `F.p.value`

field contains p-values for the F-test on *all* coefficients, including the intercept - this corresponds to the null hypothesis that all expression values are zero, and is unlikely to be sensible for single-channel arrays (and probably of limited utility for two-channel arrays).

Your motivation for using the F-test is also misguided. The F-test involving all terms is not uniformly more powerful than a t-test for a specific term, so your pre-filtering might end up throwing out genes that would otherwise be detected by using the t-test directly.

```
set.seed(42)
y <- matrix(rnorm(60000), ncol=10)
y[,1:2] <- y[,1:2] + 2 # DE genes in group 1.
library(limma)
design <- model.matrix(~gl(5,2))
fit <- lmFit(y, design)
fit <- eBayes(fit)
res.F <- topTable(fit, n=Inf, sort.by="n")
res.t <- topTable(fit, coef=2, n=Inf, sort.by="n")
# Lots of genes lost by the F-test that are
# detected directly by the t-test.
table(SigF=res.F$adj.P.Val <= 0.05,
SigT=res.t$adj.P.Val <= 0.05)
```

In addition, your second round of multiple testing correction fails to account for the fact that the gene list has effectively been pre-filtered based on the same data used in the t-test! This undermines the validity of the multiple testing, resulting in more false positives in the t-test DE gene list.

```
set.seed(1000)
y <- matrix(rnorm(60000), ncol=6)
# Adding lots of DE genes, but only in group 3.
N <- 2000
y[1:N,5:6] <- y[1:N,5:6] + 2
# Doing an F-test for any difference between groups.
library(limma)
design <- model.matrix(~gl(3,2))
fit <- lmFit(y, design)
fit <- eBayes(fit)
res.F <- topTable(fit, n=Inf, sort.by="n")
# Subsetting and repeating for the 2 v 1 comparison
keep <- res.F$adj.P.Val <= 0.05
res.2 <- topTable(fit[keep,], coef=2, n=Inf, sort.by="n")
# Consistently non-zero detected features, despite the fact that
# the null is true for all 2 v 1 comparisons - these are false positives.
sum(res.2$adj.P.Val <= 0.05)
# Correct FDR adjustment should only yield non-zero values 5% of the time.
unfilt <- topTable(fit, coef=2, n=Inf, sort.by="n")
sum(unfilt$adj.P.Val <= 0.05)
```

I would suggest *not* trying to do something fancy, and just to do the test that corresponds directly to your scientific question. If you're interested in the pairwise comparisons - just do them.