Question: GSVA enrichment score and heatmap
0
gravatar for steve.booth
9 weeks ago by
steve.booth10
steve.booth10 wrote:

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

gsea gsva R ssgsea • 250 views
ADD COMMENTlink modified 9 weeks ago by Kevin Blighe9.3k • written 9 weeks ago by steve.booth10
2
gravatar for Kevin Blighe
9 weeks ago by
Kevin Blighe9.3k
London
Kevin Blighe9.3k wrote:

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

heat
image uploader

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

gsva
free image hosting

ADD COMMENTlink modified 9 weeks ago • written 9 weeks ago by Kevin Blighe9.3k
1

Thanks so much Kevin, works perfectly for my data also.

Best, Steve

ADD REPLYlink written 9 weeks ago by steve.booth10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1394 users visited in the last hour