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

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)