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"))
Thanks so much Kevin, works perfectly for my data also.
Best, Steve
Hi Kevin, sorry to follow up after such a long time. I tried running GSVA again recently for some new data but I am getting a strange error when trying to access the es.obs using exprs(). Also when I run your examples above using the package data I receive the same error now:
Any ideas?? Thanks again, Steve
Hey Steve. Yes, between now and when I originally posted, the behaviour of the
gsva()
function changed. So, runninggsva()
now should produce the expression matrix directly and de-necessitate the need to access it via$es.obs
Ohh I see, its now under @assayData$exprs.
Thanks for the quick response! Steve