### R code from vignette source 'vignettes/GSVA/inst/doc/GSVA.Rnw' ################################################### ### code chunk number 1: options ################################################### options(width=60) ################################################### ### code chunk number 2: GSVA.Rnw:267-272 ################################################### library(GSEABase) library(GSVAdata) data(c2BroadSets) c2BroadSets ################################################### ### code chunk number 3: GSVA.Rnw:277-285 ################################################### library(Biobase) library(genefilter) library(limma) library(RColorBrewer) library(RBGL) library(graph) library(Rgraphviz) library(GSVA) ################################################### ### code chunk number 4: GSVA.Rnw:291-293 ################################################### cacheDir <- system.file("extdata", package="GSVA") cachePrefix <- "cache4vignette_" ################################################### ### code chunk number 5: GSVA.Rnw:299-300 (eval = FALSE) ################################################### ## file.remove(paste(cacheDir, list.files(cacheDir, pattern=cachePrefix), sep="/")) ################################################### ### code chunk number 6: GSVA.Rnw:322-326 ################################################### data(leukemia) leukemia_eset head(pData(leukemia_eset)) table(leukemia_eset$subtype) ################################################### ### code chunk number 7: figIQR ################################################### png(filename="GSVA-figIQR.png", width=500, height=500, res=150) IQRs <- esApply(leukemia_eset, 1, IQR) plot.ecdf(IQRs, pch=".", xlab="Interquartile range (IQR)", main="Leukemia data") abline(v=quantile(IQRs, prob=0.5), lwd=2, col="red") dev.off() ################################################### ### code chunk number 8: GSVA.Rnw:354-359 ################################################### 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") filtered_eset leukemia_filtered_eset <- filtered_eset$eset ################################################### ### code chunk number 9: GSVA.Rnw:372-375 ################################################### cache(leukemia_es <- gsva(leukemia_filtered_eset, c2BroadSets, min.sz=10, max.sz=500, verbose=TRUE)$es.obs, dir=cacheDir, prefix=cachePrefix) ################################################### ### code chunk number 10: GSVA.Rnw:389-391 ################################################### adjPvalueCutoff <- 0.001 logFCcutoff <- log2(2) ################################################### ### code chunk number 11: GSVA.Rnw:396-405 ################################################### design <- model.matrix(~ factor(leukemia_es$subtype)) colnames(design) <- c("ALL", "MLLvsALL") fit <- lmFit(leukemia_es, 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) ################################################### ### code chunk number 12: GSVA.Rnw:411-421 ################################################### logFCcutoff <- log2(2) design <- model.matrix(~ factor(leukemia_eset$subtype)) colnames(design) <- c("ALL", "MLLvsALL") fit <- lmFit(leukemia_filtered_eset, design) fit <- eBayes(fit) allGenes <- topTable(fit, coef="MLLvsALL", number=Inf) DEgenes <- topTable(fit, coef="MLLvsALL", number=Inf, p.value=adjPvalueCutoff, adjust="BH", lfc=logFCcutoff) res <- decideTests(fit, p.value=adjPvalueCutoff, lfc=logFCcutoff) summary(res) ################################################### ### code chunk number 13: leukemiaVolcano ################################################### png(filename="GSVA-leukemiaVolcano.png", width=800, height=500) par(mfrow=c(1,2)) plot(allGeneSets$logFC, -log10(allGeneSets$P.Value), pch=".", cex=4, col=grey(0.75), main="Gene sets", xlab="GSVA enrichment score difference", ylab=expression(-log[10]~~~Raw~P-value)) abline(h=-log10(max(allGeneSets$P.Value[allGeneSets$adj.P.Val <= adjPvalueCutoff])), col=grey(0.5), lwd=1, lty=2) points(allGeneSets$logFC[match(DEgeneSets$ID, allGeneSets$ID)], -log10(allGeneSets$P.Value[match(DEgeneSets$ID, allGeneSets$ID)]), pch=".", cex=4, col="red") text(max(allGeneSets$logFC)*0.85, -log10(max(allGeneSets$P.Value[allGeneSets$adj.P.Val <= adjPvalueCutoff])), sprintf("%.1f%% FDR", 100*adjPvalueCutoff), pos=1) plot(allGenes$logFC, -log10(allGenes$P.Value), pch=".", cex=4, col=grey(0.75), main="Genes", xlab="Log fold-change", ylab=expression(-log[10]~~~Raw~P-value)) abline(h=-log10(max(allGenes$P.Value[allGenes$adj.P.Val <= adjPvalueCutoff])), col=grey(0.5), lwd=1, lty=2) abline(v=c(-logFCcutoff, logFCcutoff), col=grey(0.5), lwd=1, lty=2) points(allGenes$logFC[match(DEgenes$ID, allGenes$ID)], -log10(allGenes$P.Value[match(DEgenes$ID, allGenes$ID)]), pch=".", cex=4, col="red") text(max(allGenes$logFC)*0.85, -log10(max(allGenes$P.Value[allGenes$adj.P.Val <= adjPvalueCutoff])), sprintf("%.1f%% FDR", 100*adjPvalueCutoff), pos=1) dev.off() ################################################### ### code chunk number 14: leukemiaHeatmapGeneSets ################################################### png(filename="GSVA-leukemiaHeatmapGeneSets.png", width=500, height=500) GSVAsco <- exprs(leukemia_es[DEgeneSets$ID, ]) colorLegend <- c("darkred", "darkblue") names(colorLegend) <- c("ALL", "MLL") sample.color.map <- colorLegend[pData(leukemia_es)[, "subtype"]] names(sample.color.map) <- colnames(GSVAsco) sampleClustering <- hclust(as.dist(1-cor(GSVAsco, method="spearman")), method="complete") geneSetClustering <- hclust(as.dist(1-cor(t(GSVAsco), method="pearson")), method="complete") heatmap(GSVAsco, ColSideColors=sample.color.map, xlab="samples", ylab="Gene sets and pathways", margins=c(2, 20), labRow=substr(gsub("_", " ", gsub("^KEGG_|^REACTOME_|^BIOCARTA_", "", rownames(GSVAsco))), 1, 35), labCol="", scale="row", Colv=as.dendrogram(sampleClustering), Rowv=as.dendrogram(geneSetClustering)) legend("topleft", names(colorLegend), fill=colorLegend, inset=0.01, bg="white") dev.off() ################################################### ### code chunk number 15: leukemiaHeatmapGenes ################################################### png(filename="GSVA-leukemiaHeatmapGenes.png", width=500, height=500) exps <- exprs(leukemia_eset[DEgenes$ID, ]) colorLegend <- c("darkred", "darkblue") names(colorLegend) <- c("ALL", "MLL") sample.color.map <- colorLegend[pData(leukemia_eset)[, "subtype"]] names(sample.color.map) <- colnames(exps) sampleClustering <- hclust(as.dist(1-cor(exps, method="spearman")), method="complete") geneClustering <- hclust(as.dist(1-cor(t(exps), method="pearson")), method="complete") heatmap(exps, ColSideColors=sample.color.map, xlab="samples", ylab="Genes", labRow="", labCol="", scale="row", Colv=as.dendrogram(sampleClustering), Rowv=as.dendrogram(geneClustering), margins=c(2,2)) legend("topleft", names(colorLegend), fill=colorLegend, inset=0.01, bg="white") dev.off() ################################################### ### code chunk number 16: GSVA.Rnw:529-536 ################################################### data(gbm_VerhaakEtAl) gbm_eset head(featureNames(gbm_eset)) table(gbm_eset$subtype) data(brainTxDbSets) sapply(brainTxDbSets, length) lapply(brainTxDbSets, head) ################################################### ### code chunk number 17: GSVA.Rnw:541-542 ################################################### gbm_es <- gsva(gbm_eset, brainTxDbSets, mx.diff=FALSE, verbose=FALSE)$es.obs ################################################### ### code chunk number 18: gbmSignature ################################################### png(filename="GSVA-gbmSignature.png", width=700, height=500) subtypeOrder <- c("Proneural", "Neural", "Classical", "Mesenchymal") sampleOrderBySubtype <- sort(match(gbm_es$subtype, subtypeOrder), index.return=TRUE)$ix subtypeXtable <- table(gbm_es$subtype) subtypeColorLegend <- c(Proneural="red", Neural="green", Classical="blue", Mesenchymal="orange") geneSetOrder <- c("astroglia_up", "astrocytic_up", "neuronal_up", "oligodendrocytic_up") geneSetLabels <- gsub("_", " ", geneSetOrder) hmcol <- colorRampPalette(brewer.pal(10, "RdBu"))(256) hmcol <- hmcol[length(hmcol):1] heatmap(exprs(gbm_es)[geneSetOrder, sampleOrderBySubtype], Rowv=NA, Colv=NA, scale="row", margins=c(3,5), col=hmcol, ColSideColors=rep(subtypeColorLegend[subtypeOrder], times=subtypeXtable[subtypeOrder]), labCol="", gbm_es$subtype[sampleOrderBySubtype], labRow=paste(toupper(substring(geneSetLabels, 1,1)), substring(geneSetLabels, 2), sep=""), cexRow=2, main=" \n ") par(xpd=TRUE) text(0.22,1.11, "Proneural", col="red", cex=1.2) text(0.36,1.11, "Neural", col="green", cex=1.2) text(0.48,1.11, "Classical", col="blue", cex=1.2) text(0.66,1.11, "Mesenchymal", col="orange", cex=1.2) mtext("Gene sets", side=4, line=0, cex=1.5) mtext("Samples ", side=1, line=4, cex=1.5) dev.off() ################################################### ### code chunk number 19: GSVA.Rnw:605-609 ################################################### canonicalC2BroadSets <- c2BroadSets[c(grep("^KEGG", names(c2BroadSets)), grep("^REACTOME", names(c2BroadSets)), grep("^BIOCARTA", names(c2BroadSets)))] canonicalC2BroadSets ################################################### ### code chunk number 20: GSVA.Rnw:617-620 ################################################### cache(leukemia_canonicalPwy_es <- gsva(leukemia_eset, canonicalC2BroadSets, min.sz=10, max.sz=500, mx.diff=TRUE, verbose=TRUE)$es.obs, dir=cacheDir, prefix=cachePrefix) ################################################### ### code chunk number 21: GSVA.Rnw:628-630 ################################################### overlapMatrix <- computeGeneSetsOverlap(canonicalC2BroadSets, leukemia_eset, min.sz=10, max.sz=500) ################################################### ### code chunk number 22: GSVA.Rnw:651-652 ################################################### library(qpgraph) ################################################### ### code chunk number 23: GSVA.Rnw:666-670 ################################################### cache(gennrr <- qpGenNrr(leukemia_canonicalPwy_es, datasetIdx="subtype", qOrders=c(ALL=15, MLL=12), identicalQs=FALSE, clusterSize=2)$genNrr, dir=cacheDir, prefix=cachePrefix) ################################################### ### code chunk number 24: GSVA.Rnw:675-676 ################################################### gennrr[overlapMatrix > 0.05] <- NA ################################################### ### code chunk number 25: gennrrmaxbd ################################################### png(filename="GSVA-gennrrmaxbd.png", width=800, height=800, res=150) qpbd <- qpBoundary(gennrr, N=ncol(leukemia_canonicalPwy_es), breaks=20) dev.off() ################################################### ### code chunk number 26: GSVA.Rnw:708-711 ################################################### g <- qpGraph(gennrr, threshold=qpbd$threshold, return.type="graphNEL") g max(degree(g)) ################################################### ### code chunk number 27: GSVA.Rnw:730-732 ################################################### pac <- qpPAC(leukemia_canonicalPwy_es, g, return.K=TRUE, tol=0.01, verbose=FALSE) ################################################### ### code chunk number 28: GSVA.Rnw:743-756 ################################################### gAM <- qpGraph(gennrr, threshold=qpbd$threshold) ridx <- row(pac$P)[as.matrix(upper.tri(pac$P) & gAM)] cidx <- col(pac$P)[as.matrix(upper.tri(pac$P) & gAM)] allEdges <- data.frame(PWYi=colnames(pac$P)[ridx], PWYj=colnames(pac$P)[cidx], PAC=pac$R[cbind(ridx, cidx)], PAC.P.value=pac$P[cbind(ridx, cidx)], PAC.adj.P.value=p.adjust(pac$P[cbind(ridx, cidx)], method="fdr"), PCC=cov2cor(solve(pac$K))[cbind(ridx, cidx)]) allEdges <- allEdges[sort(allEdges$PAC.adj.P.value, index.return=TRUE)$ix, ] dim(allEdges) ################################################### ### code chunk number 29: GSVA.Rnw:762-764 ################################################### sigEdges <- allEdges[which(allEdges$PAC.adj.P.value < 0.1 & abs(allEdges$PCC) > 0.7) , ] dim(sigEdges) ################################################### ### code chunk number 30: GSVA.Rnw:769-775 ################################################### vtc <- unique(as.character(unlist(sigEdges[, c("PWYi", "PWYj")], use.names=FALSE))) g <- new("graphNEL", nodes=vtc, edgemode="undirected") g <- addEdge(from=as.character(sigEdges[, "PWYi"]), to=as.character(sigEdges[, "PWYj"]), graph=g) g ################################################### ### code chunk number 31: leukemiaNet ################################################### png(filename="GSVA-leukemiaNet.png", width=800, height=800, res=150) nodlab <- gsub("_", " ", gsub("KEGG_|REACTOME_|BIOCARTA_", "", vtc)) nodlab <- sapply(nodlab, function(x) { v <- unlist(strsplit(x, ' ')) ; t <- ""; l <- 0; for (w in v) { t <- paste(t,w,sep=" ") ; if (nchar(t)-l > 5) { t <- paste(t, "\ \n", sep="") ; l <- nchar(t) } } ; t }) names(nodlab) <- vtc nodeRenderInfo(g) <- list(shape="ellipse", label=nodlab, fill="lightgrey", lwd=1) edgeRenderInfo(g) <- list(lwd=1) g <- layoutGraph(g, layoutType="fdp") renderGraph(g) dev.off() ################################################### ### code chunk number 32: GSVA.Rnw:800-802 ################################################### cct <- table(listLen(connComp(g))) cct ################################################### ### code chunk number 33: GSVA.Rnw:810-814 ################################################### vtcCCmax <- connComp(g)[[which(sapply(connComp(g), length) == max(as.integer(names(cct))))]] gCCmax <- subGraph(vtcCCmax, g) gCCmax table(degree(gCCmax)) ################################################### ### code chunk number 34: GSVA.Rnw:820-821 ################################################### sort(degree(gCCmax), decreasing=TRUE)[1:4] ################################################### ### code chunk number 35: GSVA.Rnw:827-836 ################################################### vtcHubNet <- names(sort(degree(gCCmax), decreasing=TRUE)[1:4]) pairs <- combn(vtcHubNet, 2) sp <- sp.between(gCCmax, pairs[1, ], pairs[2, ]) vtcInShortestPaths <- unlist(sapply(sp, function(x) x$path_detail), use.names=FALSE) vtcHubNet <- unique(c(vtcHubNet, unlist(boundary(vtcHubNet, gCCmax), use.names=FALSE), vtcInShortestPaths)) gHubNet <- subGraph(vtcHubNet, gCCmax) gHubNet ################################################### ### code chunk number 36: selectedGGMhighDeg ################################################### nodlab <- gsub("_", " ", gsub("KEGG_|REACTOME_|BIOCARTA_", "", vtcHubNet)) nodlab <- sapply(nodlab, function(x) { v <- unlist(strsplit(x, ' ')) ; t <- ""; l <- 0; for (w in v) { t <- paste(t,w,sep=" ") ; if (nchar(t)-l > 5) { t <- paste(t, "\ \n", sep="") ; l <- nchar(t) } } ; t }) names(nodlab) <- vtcHubNet nodeRenderInfo(gHubNet) <- list(shape="ellipse", label=nodlab, fill="lightgrey", lwd=1) edgeRenderInfo(gHubNet) <- list(lwd=1) gHubNet <- layoutGraph(gHubNet, layoutType="fdp") renderGraph(gHubNet) ################################################### ### code chunk number 37: GSVA.Rnw:881-885 ################################################### ids <- geneIds(c2BroadSets["KEGG_ONE_CARBON_POOL_BY_FOLATE"])[[1]] unlist(mget(ids[!is.na(match(ids, unlist(mget(featureNames(leukemia_eset), hgu95aENTREZID))))], org.Hs.egSYMBOL), use.names=FALSE) ################################################### ### code chunk number 38: info ################################################### toLatex(sessionInfo())