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
**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()

<a href="https://imgbb.com/"><img src="https://i.ibb.co/tJFMrMG/heat.png" alt="heat" border="0"></a>

#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"))

<a href="https://ibb.co/Q8DZdFD"><img src="https://i.ibb.co/1sz36fz/gsva.png" alt="gsva" border="0"></a>


