Question: GSVA enrichment score and heatmap
2
gravatar for steve.booth
18 months ago by
steve.booth30
steve.booth30 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 • 2.3k views
ADD COMMENTlink modified 18 months ago by Kevin Blighe41k • written 18 months ago by steve.booth30
2
gravatar for Kevin Blighe
18 months ago by
Kevin Blighe41k
Guy's Hospital, London
Kevin Blighe41k 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 18 months ago • written 18 months ago by Kevin Blighe41k
1

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

Best, Steve

ADD REPLYlink written 18 months ago by steve.booth30
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: 1638 users visited in the last hour