Dear all,

I am trying to use the limma ROAST and camera() to conduct enrichment analysis. However, the two method gave me different results, and the ROAST results were pretty strange. It showed that all pathways were upregulated. Even when I changed another expression matrix, the Pvalue was still the same. I am wondering where I was wrong.

Here are my scripts.

expr.i <- exprCPM %>%
  select(contains("0h") | contains("0.5h"))
target.i <- targets[colnames(expr.i),]
design <- model.matrix(~ 0 + Time, data = target.i)

#voom <- voom(expr.i, design, plot = TRUE)

index1 <- rownames(expr.i) %in% AP1
index2 <- rownames(expr.i) %in% NRF2
index3 <- rownames(expr.i) %in% UPR
index4 <- rownames(expr.i) %in% HIF1
index5 <- rownames(expr.i) %in% NFkappaB

GeneList <- list(AP1  = index1,
                 NRF2 = index2,
                 UPR  = index3,
                 HIF1 = index4,
                 NFkappaB = index5)

camera(expr.i, GeneList, design = design, contrast = 2)
#         NGenes Direction      PValue        FDR
#AP1          45        Up 0.003845355 0.01059174
#UPR          26        Up 0.004236695 0.01059174
#NRF2         22        Up 0.035639490 0.05939915
#HIF1         34        Up 0.084318148 0.10539768
#NFkappaB     33        Up 0.344350372 0.34435037

mroast(expr.i, GeneList, design = design, contrast = 2)
#         NGenes   PropDown    PropUp Direction PValue   FDR PValue.Mixed FDR.Mixed
#UPR          26 0.03846154 0.9615385        Up  5e-04 5e-04        5e-04     5e-04
##AP1          45 0.08888889 0.9111111        Up  5e-04 5e-04        5e-04     5e-04
#NRF2         22 0.09090909 0.9090909        Up  5e-04 5e-04        5e-04     5e-04
#HIF1         34 0.14705882 0.8529412        Up  5e-04 5e-04        5e-04     5e-04
#NFkappaB     33 0.18181818 0.8181818        Up  5e-04 5e-04        5e-04     5e-04

Here is my data used in ROAST.

 > head(expr.i)
      EGI1_0h_1 EGI1_0h_2 EGI1_0h_3 EGI1_0.5h_1 EGI1_0.5h_2 EGI1_0.5h_3
7105   4.477632  4.702390  4.941865    4.673160    4.763974    4.857584
8813   5.386476  5.427719  5.606738    5.434388    5.455882    5.581303
57147  3.285559  3.048268  3.185166    3.040048    3.117499    3.157894
55732  5.160325  5.217473  5.212180    5.159004    5.262421    5.280877
2268  -2.224361 -2.550945 -2.550945   -2.550945   -2.550945   -2.550945
3075   2.115124  1.884813  2.070823    2.214104    2.054585    1.813567

> design
            Time0.5h Time0h
EGI1_0h_1          0      1
EGI1_0h_2          0      1
EGI1_0h_3          0      1
EGI1_0.5h_1        1      0
EGI1_0.5h_2        1      0
EGI1_0.5h_3        1      0
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$Time
[1] "contr.treatment"

Thanks for answer! All the best.

The contrast = 2 in an intercept-less design test whether the expression in the 2nd group is not 0. You can easily solve this by using ~time as design given that you only have two groups here. Alternatively, use the output of makeContrasts().