#### The support.bioconductor.org editor has been updated to markdown! Please see more info at: Tutorial: Updated Support Site Editor

Question: limma roast syntax and interpretation
1
4.3 years ago by
Juliet Hannah360
United States
Juliet Hannah360 wrote:

I am using ROAST. I am a little confused by the results, and I wanted to check if my syntax appears correct. I understand that geneSetTest is using a different method than ROAST.  However, I am obtaining substantially different results. Also, I am how does one interpret "Active.Prop=1". Thanks.

# limma syntax

# stage has four levels, and I am interested in the difference between B and C

design <- model.matrix(~0+key$stage) colnames(design) = c("A","B","C","D") head(design) # A 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 fit <- lmFit(e, design) fit <- eBayes(fit) cont.matrix <- makeContrasts(risk="B-C",levels=design) fit2 <- contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2) tt2 <-topTable(fit2, number=Inf,coef=1) tt2$PROBEID = rownames(tt2)

# use geneSetTest

# read in probes (not shown)
# tt2 is the topTable

myIndex = rownames(tt2) %in% probes
geneSetTest(myIndex, tt2$t) # 0.2439185 # use roast myIndex = rownames(e) %in% probes roast(e,myIndex,design,contrasts=cont.matrix,nrot=5000) # the result differs substantially from geneSetTest # what does it mean that "Active.Prop=1" # Active.Prop P.Value #Down 0 1.00000000 #Up 1 0.00019996 #Mixed 1 0.00019996 limma roast • 2.0k views ADD COMMENTlink modified 4.3 years ago by Gordon Smyth36k • written 4.3 years ago by Juliet Hannah360 Answer: limma roast syntax and interpretation 2 4.3 years ago by Gordon Smyth36k Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia Gordon Smyth36k wrote: geneSetTest and roast test different hypotheses. As the help pages says: "geneSetTest performs a competitive test in the sense that genes in the test set are compared to other genes. If the statistic is a genewise test statistic for differential expression, then geneSetTest tests whether genes in the set are more differentially expressed than genes not in the set. By contrast, a self-contained gene set test such as roast tests whether genes in the test set are differentially expressed, in an absolute sense, without regard to any other genes on the array." In your case, the roast test tells you that your gene set is up-regulated. In fact every single gene in the set has a logFC more than one standard error above zero (hence the active proportion is 100%). The geneSetTest test tells you however than this gene set does not stand out from the other genes. It would seem that most or all of the genes in the fit have positive t-statistics about the same size as those for genes in your specified set. This does beg the question as to which genes are included in your matrix e. geneSetTest only makes sense if e contains all expressed genes. It is only possible for all genes to be up in B vs C if the data is not normalized properly or you have specially selected genes that are DE. BTW, it would be clearer if you defined myIndex just once rather than in two different ways. Then you would be sure that the two tests are working with the same set: cont.matrix <- makeContrasts(risk="B-C",levels=design) fit2 <- contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2) myIndex = rownames(e) %in% probes geneSetTest(myIndex, fit2$t, alternative="down")
geneSetTest(myIndex, fit2\$t, alternative="up")
roast(e, myIndex, design, contrast=cont.matrix, nrot=5000)

Thanks, as always, for your help. In this post ( limma roast syntax for overall anova ), I had a question about ROAST, and your answer directed me to geneSetTest, which I misinterpreted to mean that one could be substituted for another. I should have read the documentation more carefully. What I did not mention in my question is that looking at the upregulated genes, there does not appear to be much of a signal, which is why I couldn't see why ROAST would indicate one. I assumed  I made a mistake somewhere.

ADD REPLYlink modified 4.3 years ago by Gordon Smyth36k • written 4.3 years ago by Juliet Hannah360

For the last line of your code:

roast(e, myIndex, design, contrasts=cont.matrix, nrot=5000)

Shouldn't it be "contrast" instead of "contrasts"? And shouldn't the argument be a numerical vector corresponding to the columns of the design matrix instead of a full contrast matrix?

1
1. Yes.
2. contrast is fine as a single-column matrix, which is taken to be a numeric vector. Equivalently you could specify contrast=cont.matrix[,1].

what exactly is e here (from "roast(e, myIndex, design, contrast=cont.matrix, nrot=5000)")?

e is the same in my answer as in the original question. It is the expression object.