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
gravatar for Juliet Hannah
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
gravatar for Gordon Smyth
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)

 

ADD COMMENTlink modified 6 months ago • written 4.3 years ago by Gordon Smyth36k

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?

ADD REPLYlink written 6 months ago by siajunren0
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].
ADD REPLYlink modified 6 months ago • written 6 months ago by Gordon Smyth36k

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

ADD REPLYlink written 27 days ago by Daniel10

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

ADD REPLYlink written 26 days ago by Gordon Smyth36k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 227 users visited in the last hour