Question: Hi there,
I have just used GSVA to perform ssGSEA + LIMMA as follows:

eDat_es <- gsva(eset, c2BroadSets, min.sz=10, max.sz=500, method="gsva", mx.diff=TRUE,verbose=FALSE)$es.obs

design <- model.matrix(~ Macrophages.M1 + Age + Sex+ packyears +BMI, data = pDat_MR)

fit <- lmFit(eDat_es, design)

fit <- eBayes(fit)

allGeneSets <- topTable(fit, coef=2, number=Inf)

M1_DEgeneSets <- topTable(fit, coef=2, number=Inf,
p.value=adjPvalueCutoff, adjust="BH")

res <- decideTests(fit, p.value=adjPvalueCutoff)

summary(res)

Does anyone know how I can access the matrix of enrichment scores that is produced by the gsva command and plot as a heatmap? This is shown in figure 4 of the package vignette but not code provided.

Thanks!
Steve

Answer: **Edit: I hadn't noticed your reference to Figure 4. I initially answered for Figure 6 but have added Figure 4**

Hey Steve, I presume that you are referring to Figure 6 from the manual? - https://www.bioconductor.org/packages/devel/bioc/vignettes/GSVA/inst/doc/GSVA.pdf

The key is the *es.obs* child object of the *gsva* parent object. *es.obs* contains the numerical values, which can be accessed with `exprs()`.

I have managed to roughly reproduce the figure from the manual as follows (for both, you'll have to add your own colour bar labels with the ColSideColors parameter of `heatmap.2`):

#Figure 4

filtered_eset <- nsFilter(leukemia_eset, require.entrez=TRUE, remove.dupEntrez=TRUE, var.func=IQR, var.filter=TRUE, var.cutoff=0.5, filterByQuantile=TRUE, feature.exclude="^AFFX")

leukemia_filtered_eset <- filtered_eset$eset

require("limma")
adjPvalueCutoff <- 0.001
logFCcutoff <- log2(4)

design <- model.matrix(~factor(leukemia_eset$subtype))
colnames(design) <- c("ALL", "MLLvsALL")
fit <- lmFit(leukemia_eset, design)
fit <- eBayes(fit)
allGeneSets <- topTable(fit, coef="MLLvsALL", number=Inf)
DEgeneSets <- topTable(fit, coef="MLLvsALL", number=Inf, p.value=adjPvalueCutoff, adjust="BH")
res <- decideTests(fit, p.value=adjPvalueCutoff)
summary(res)

sig.genes <- c(names(res[res[,2]==1,2]), names(res[res[,2]==-1,2]))
leukemia_eset.siggenes <- leukemia_eset[which(rownames(leukemia_eset) %in% sig.genes),]

leukemia_es <- gsva(leukemia_eset.siggenes, c2BroadSets, min.sz=10, max.sz=500, verbose=TRUE)

pdf("heat.pdf", width=5, height=8)
heatmap.2(exprs(leukemia_es$es.obs), dendrogram="col", main="", key=FALSE, scale="none", density.info="none", reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean), trace="none", cexRow=0.2, labCol=NA, distfun=function(x) as.dist(1-cor(t(x))), hclustfun=function(x) hclust(x, method="ward.D2"))
dev.off()

#Figure 6
gbm_es <- gsva(gbm_eset, brainTxDbSets, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)

heat <- exprs(gbm_es$es.obs)
heat <- heat[grep("_up", rownames(heat)), ]
heat <- t(scale(t(heat)))

require(rColorBrewer)
myCol <- colorRampPalette(c("darkblue", "white", "red4"))(100)
myBreaks <- seq(-2, 2, length.out=101)

require("gplots")

rownames(heat) <- gsub("_up", "\nup", rownames(heat))

par(mar=c(2,2,2,2), cex=0.8)
heatmap.2(heat, col=myCol, breaks=myBreaks, dendrogram="none", main="", key=FALSE, scale="none", density.info="none", reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean), trace="none", cexRow=0.65, labCol=NA, distfun=function(x) as.dist(1-cor(t(x))), hclustfun=function(x) hclust(x, method="ward.D2"))