Question: GSVA enrichment score and heatmap
2
gravatar for steve.booth
2.4 years ago by
steve.booth40
steve.booth40 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 • 3.9k views
ADD COMMENTlink modified 2.4 years ago by Kevin Blighe55k • written 2.4 years ago by steve.booth40
2
gravatar for Kevin Blighe
2.4 years ago by
Kevin Blighe55k
Kevin Blighe55k 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 2.4 years ago • written 2.4 years ago by Kevin Blighe55k
1

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

Best, Steve

ADD REPLYlink written 2.4 years ago by steve.booth40

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:

    data(gbm_VerhaakEtAl)
    data(brainTxDbSets)
    gbm_es <- gsva(gbm_eset, brainTxDbSets, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)  
    heat <- exprs(gbm_es$es.obs)
Error in (function (classes, fdef, mtable)  : 
      unable to find an inherited method for function ‘exprs’ for signature ‘"NULL"’

Any ideas?? Thanks again, Steve

ADD REPLYlink modified 9 months ago • written 9 months ago by steve.booth40
1

Hey Steve. Yes, between now and when I originally posted, the behaviour of the gsva() function changed. So, running gsva() now should produce the expression matrix directly and de-necessitate the need to access it via $es.obs

ADD REPLYlink modified 9 months ago • written 9 months ago by Kevin Blighe55k

Ohh I see, its now under @assayData$exprs.

Thanks for the quick response! Steve

ADD REPLYlink written 9 months ago by steve.booth40
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: 1107 users visited in the last hour