limma roast syntax and interpretation
1
1
Entering edit mode
Juliet Hannah ▴ 360
@juliet-hannah-4531
Last seen 4.9 years ago
United States

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 • 3.8k views
ADD COMMENT
2
Entering edit mode
@gordon-smyth
Last seen 10 hours ago
WEHI, Melbourne, Australia

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 COMMENT
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
1
Entering edit mode
  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 REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 891 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6