limma roast syntax for overall anova
1
0
Entering edit mode
Juliet Hannah ▴ 360
@juliet-hannah-4531
Last seen 3.2 years ago
United States

All,

I am testing an overall anova (4 groups; is there any difference across groups) using the following design matrix in limma

  (Intercept) B C D
1           1 0 0 0
2           1 0 0 0
3           1 0 0 0
4           1 0 0 0
5           1 0 0 0
6           1 0 0 0

I have used ROAST (in limma) in the past with a specific contrast, such as

roast(expression_matrix,index,design_matrix,contrast=10,nrot=3000).

In the current case, I didn't need any contrasts. In limma, I had used

design <- model.matrix(~key$stage) colnames(design)[2:4] = c("B","C","D") fit <- lmFit(e, design) fit <- eBayes(fit) topTable(fit,adjust.method="BH") to obtain the overall F-test (I understand this test is not that interesting, but it is requested of me). Using this setup, how can I specify this overall test within ROAST. Thanks for your help. Juliet limma • 1.9k views ADD COMMENT 3 Entering edit mode @gordon-smyth Last seen 1 minute ago WEHI, Melbourne, Australia Dear Juliet, To get the overall anova F-test you need:  topTable(fit[,-1]) or alternatively  topTable(fit, coef=2:4) which gives the same result. Note that the intercept needs to be removed from the test, as it always is in classical anova, and that's what the -1 does. You also ask "Using this setup, how can I specify this overall test within ROAST", and that is a good question. We haven't implemented F-tests for ROAST. If we get enough requests, we will consider doing so, but we haven't so far. If the meantime, you can input the F-statistics to the geneSetTest() function:  fit2 <- fit[,-1] geneSetTest(index, fit2$F)

etc.

Best wishes
Gordon

0
Entering edit mode

Hi Gordon,

Thanks for your help (and software!).

I confusingly showed only the first few rows of the design matrix. I was trying to convey that I left the intercept in.

limma gave me the message:

> mytt <- topTable(fit,number=Inf)
Removing intercept from test coefficients

which gives me the same result as the more clear version you gave me

 mytt <- topTable(fit,coef=2:4,number=Inf).

I will proceed with geneSetTest.

Thanks,

Juliet

0
Entering edit mode

Oh, yes, that's right. I had forgotten that limma now recognizes the special intercept column name "(Intercept)" and removes it automatically. So you already had the correct anova test.

Best wishes
Gordon