Hi,

I'm analysing some microarray and rnaseq data for gene set enrichment analysis. I have three independent sample groups. I'm using GSVA to do the enrichment analysis then using limma to fit a model to the GSVA enrichment scores.

library(GSVA)

library(limma)

p <- 10 ## number of genes

n <- 45 ## number of samples

nGrp1 <- 15 ## number of samples in group 1

nGrp2 <- 15 ## number of samples in group 2

nGrp3 <- 15 ## number of samples in group 2

## consider three disjoint gene sets

geneSets <- list(set1=paste("g", 1:3, sep=""),

set2=paste("g", 4:6, sep=""),

set3=paste("g", 7:10, sep=""))

## sample data from a normal distribution with mean 0 and st.dev. 1

y <- matrix(rnorm(n*p), nrow=p, ncol=n,

dimnames=list(paste("g", 1:p, sep="") , paste("s", 1:n, sep="")))

## genes in set1 are expressed at higher levels in the last 10 samples

y[geneSets$set1, (nGrp1+1):n] <- y[geneSets$set1, (nGrp1+1):n] + 2

## build design matrix

design <- cbind(intercept=1, groups=c(rep(0, nGrp1), rep(1, nGrp2), rep(2, nGrp3)))

## estimate GSVA enrichment scores for the three sets

gsva_es <- gsva(y, geneSets, mx.diff=1)$es.obs

## fit linear model to the GSVA enrichment scores

fit <- lmFit(gsva_es, design)

## estimate moderated t-statistics

fit <- eBayes(fit)

## set1 is differentially expressed

topTable(fit, coef="groups")

This will give me the test-wise enrichment values and p-value, however I would like to do all additional group comparisons, e.g., between group 1vs2, 1vs3, and 2vs3. Is there a way to go about doing this with limma?

Thanks for the help.

Martin