###################################################
### chunk number 1: options
###################################################
#line 42 "vignettes/GSVA/inst/doc/GSVA.Rnw"
options(width=60)


###################################################
### chunk number 2: 
###################################################
#line 255 "vignettes/GSVA/inst/doc/GSVA.Rnw"
library(GSEABase)
library(GSVAdata)

data(c2BroadSets)
c2BroadSets


###################################################
### chunk number 3: 
###################################################
#line 265 "vignettes/GSVA/inst/doc/GSVA.Rnw"
library(Biobase)
library(genefilter)
library(limma)
library(RColorBrewer)
library(graph)
library(Rgraphviz)
library(GSVA)


###################################################
### chunk number 4: 
###################################################
#line 278 "vignettes/GSVA/inst/doc/GSVA.Rnw"
cacheDir <- system.file("extdata", package="GSVA")
cachePrefix <- "cache4vignette_"


###################################################
### chunk number 5:  eval=FALSE
###################################################
## #line 286 "vignettes/GSVA/inst/doc/GSVA.Rnw"
## file.remove(paste(cacheDir, list.files(cacheDir, pattern=cachePrefix), sep="/"))


###################################################
### chunk number 6: 
###################################################
#line 309 "vignettes/GSVA/inst/doc/GSVA.Rnw"
data(leukemia)
leukemia_eset
head(pData(leukemia_eset))
table(leukemia_eset$subtype)


###################################################
### chunk number 7: figIQR
###################################################
#line 321 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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()


###################################################
### chunk number 8: 
###################################################
#line 341 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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


###################################################
### chunk number 9: 
###################################################
#line 359 "vignettes/GSVA/inst/doc/GSVA.Rnw"
cache(leukemia_es <- gsva(leukemia_filtered_eset, c2BroadSets,
                           min.sz=10, max.sz=500, verbose=FALSE)$es.obs,
                           dir=cacheDir, prefix=cachePrefix)


###################################################
### chunk number 10: 
###################################################
#line 376 "vignettes/GSVA/inst/doc/GSVA.Rnw"
adjPvalueCutoff <- 0.001
logFCcutoff <- log2(2)


###################################################
### chunk number 11: 
###################################################
#line 383 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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)


###################################################
### chunk number 12: 
###################################################
#line 398 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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)


###################################################
### chunk number 13: leukemiaVolcano
###################################################
#line 415 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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()


###################################################
### chunk number 14: leukemiaHeatmapGeneSets
###################################################
#line 458 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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()


###################################################
### chunk number 15: leukemiaHeatmapGenes
###################################################
#line 483 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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()


###################################################
### chunk number 16: 
###################################################
#line 519 "vignettes/GSVA/inst/doc/GSVA.Rnw"
data(gbm_VerhaakEtAl)
gbm_eset
head(featureNames(gbm_eset))
table(gbm_eset$subtype)
data(brainTxDbSets)
sapply(brainTxDbSets, length)
lapply(brainTxDbSets, head)


###################################################
### chunk number 17: 
###################################################
#line 531 "vignettes/GSVA/inst/doc/GSVA.Rnw"
gbm_es <- gsva(gbm_eset, brainTxDbSets, mx.diff=FALSE, verbose=FALSE)$es.obs


###################################################
### chunk number 18: gbmSignature
###################################################
#line 546 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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()


###################################################
### chunk number 19: 
###################################################
#line 594 "vignettes/GSVA/inst/doc/GSVA.Rnw"
KEGGc2BroadSets <- c2BroadSets[grep("^KEGG", names(c2BroadSets))]
KEGGc2BroadSets


###################################################
### chunk number 20: 
###################################################
#line 604 "vignettes/GSVA/inst/doc/GSVA.Rnw"
leukemiaKEGG_es <- gsva(leukemia_eset, KEGGc2BroadSets, min.sz=10, max.sz=500,
                        mx.diff=TRUE, verbose=FALSE)$es.obs


###################################################
### chunk number 21: 
###################################################
#line 614 "vignettes/GSVA/inst/doc/GSVA.Rnw"
overlapMatrix <- computeGeneSetsOverlap(KEGGc2BroadSets, leukemia_eset,
                                        min.sz=10, max.sz=500)


###################################################
### chunk number 22: 
###################################################
#line 625 "vignettes/GSVA/inst/doc/GSVA.Rnw"
pcc <- cor(t(exprs(leukemiaKEGG_es)))
pcc[overlapMatrix > .05] <- 0
pcc[lower.tri(pcc)] <- 0
diag(pcc) <- 0
arrIdxs <- which(abs(pcc) > 0.8, arr.ind=TRUE)

pccEdges <- data.frame(PWYi=featureNames(leukemiaKEGG_es)[arrIdxs[,1]],
                       PWYj=featureNames(leukemiaKEGG_es)[arrIdxs[,2]],
                       PCC=pcc[arrIdxs])


###################################################
### chunk number 23: leukemiaPCCnet
###################################################
#line 637 "vignettes/GSVA/inst/doc/GSVA.Rnw"
png(filename="GSVA-leukemiaPCCnet.png", width=1100, height=1100, res=150)
vtc <- unique(as.character(unlist(pccEdges[, c("PWYi", "PWYj")], use.names=FALSE)))
g2plot <- new("graphNEL", nodes=vtc, edgemode="undirected")
g2plot <- addEdge(from=as.character(pccEdges[, "PWYi"]),
                  to=as.character(pccEdges[, "PWYj"]),
                  graph=g2plot)
nodlab <- gsub("_", " ", gsub("KEGG_", "", 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(g2plot) <- list(shape="ellipse", label=nodlab, fill="lightgrey", lwd=2)
g2plot <- layoutGraph(g2plot, layoutType="twopi")
renderGraph(g2plot)
dev.off()


###################################################
### chunk number 24: 
###################################################
#line 671 "vignettes/GSVA/inst/doc/GSVA.Rnw"
library(qpgraph)


###################################################
### chunk number 25: 
###################################################
#line 682 "vignettes/GSVA/inst/doc/GSVA.Rnw"
cache(avgnrr <- qpAvgNrr(leukemiaKEGG_es, verbose=FALSE),
                         dir=cacheDir, prefix=cachePrefix)


###################################################
### chunk number 26: 
###################################################
#line 689 "vignettes/GSVA/inst/doc/GSVA.Rnw"
avgnrr[overlapMatrix > 0.05] <- NA


###################################################
### chunk number 27: avgnrrcliquenr
###################################################
#line 702 "vignettes/GSVA/inst/doc/GSVA.Rnw"
png(filename="GSVA-avgnrrcliquenr.png", width=800, height=800, res=150)


###################################################
### chunk number 28: avgnrrcliquenr
###################################################
#line 705 "vignettes/GSVA/inst/doc/GSVA.Rnw"
qpclq <- qpClique(avgnrr, N=dim(leukemiaKEGG_es)[2])


###################################################
### chunk number 29: avgnrrcliquenr
###################################################
#line 708 "vignettes/GSVA/inst/doc/GSVA.Rnw"
dev.off()


###################################################
### chunk number 30: 
###################################################
#line 726 "vignettes/GSVA/inst/doc/GSVA.Rnw"
g <- qpGraph(avgnrr, threshold=qpclq$threshold)
w <- qpCliqueNumber(g, verbose=FALSE)
w


###################################################
### chunk number 31: 
###################################################
#line 748 "vignettes/GSVA/inst/doc/GSVA.Rnw"
cache(pac <- qpPAC(leukemiaKEGG_es, g, return.K=TRUE, verbose=FALSE),
                   dir=cacheDir, prefix=cachePrefix)


###################################################
### chunk number 32: 
###################################################
#line 760 "vignettes/GSVA/inst/doc/GSVA.Rnw"
ridx <- row(pac$P)[as.matrix(upper.tri(pac$P) & g)]
cidx <- col(pac$P)[as.matrix(upper.tri(pac$P) & g)]
sigEdges <- which(p.adjust(pac$P[cbind(ridx, cidx)], method="holm") < 0.05)

sigEdges <- data.frame(PWYi=colnames(pac$P)[ridx][sigEdges],
                       PWYj=colnames(pac$P)[cidx][sigEdges],
                       PAC=pac$R[cbind(ridx, cidx)][sigEdges],
                       P.value=pac$P[cbind(ridx, cidx)][sigEdges],
                       PCC=cov2cor(solve(pac$K))[cbind(ridx, cidx)][sigEdges])
sigEdges <- sigEdges[sort(abs(sigEdges$P.value), index.return=TRUE)$ix, ]
dim(sigEdges)


###################################################
### chunk number 33: leukemiaNet
###################################################
#line 774 "vignettes/GSVA/inst/doc/GSVA.Rnw"
png(filename="GSVA-leukemiaNet.png", width=800, height=800, res=150)
vtc <- unique(as.character(unlist(sigEdges[, c("PWYi", "PWYj")], use.names=FALSE)))
g2plot <- new("graphNEL", nodes=vtc, edgemode="undirected")
g2plot <- addEdge(from=as.character(sigEdges[, "PWYi"]),
                  to=as.character(sigEdges[, "PWYj"]),
                  graph=g2plot)
nodlab <- gsub("_", " ", gsub("KEGG_", "", 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(g2plot) <- list(shape="ellipse", label=nodlab, fill="lightgrey", lwd=2)
g2plot <- layoutGraph(g2plot, layoutType="neato")
renderGraph(g2plot)
dev.off()


###################################################
### chunk number 34: 
###################################################
#line 806 "vignettes/GSVA/inst/doc/GSVA.Rnw"
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)


###################################################
### chunk number 35: info
###################################################
#line 818 "vignettes/GSVA/inst/doc/GSVA.Rnw"
toLatex(sessionInfo())