DEGraph.R
In DEGraph: Two-sample tests on a graph

### R code from vignette source 'DEGraph.Rnw'

###################################################
### code chunk number 1: lib
###################################################
library(DEGraph)


###################################################
### code chunk number 2: morelib
###################################################
library("R.utils")
##library(graph)
##library(rrcov) ## for 'T2.test'
library(corpcor)
library(KEGGgraph)
library(Rgraphviz)
##library(RBGL)
library(fields) # For image.plot called in plotValuedGraph
library(lattice)
library(marray)
verbose <- TRUE


###################################################
### code chunk number 3: data
###################################################
data("Loi2008_DEGraphVignette", package="DEGraph")
classData <- classLoi2008
exprData <- exprLoi2008
annData <- annLoi2008
grList <- grListKEGG


###################################################
### code chunk number 4: ttests
###################################################
## Individual t-test p-values
X1 <- t(exprData[, classData==0])
X2 <- t(exprData[, classData==1])
ttpv <- c()
tts <- c()
for(i in 1:ncol(X1)) {
  tt <- t.test(X1[,i],X2[,i])
  ttpv[i]=unlist(tt$p.value)
  tts[i]=unlist(tt$statistic)
}
names(ttpv) <- names(tts) <- rownames(exprData)


###################################################
### code chunk number 5: mvtests
###################################################
prop <- 0.2
## Multivariate tests
resList <- NULL
for (ii in seq(along=grList)) {
  gr <- grList[[ii]]
  res <- testOneGraph(gr, exprData, classData, verbose=verbose, prop=prop)
  resList <- c(resList, list(res))
}
resNames <- names(grList)
pLabels <- sapply(grList, attr, "label")

## get rid of NULL results (no connected component of size > 1)
isNULL <- sapply(resList, is.null)
if (sum(isNULL)) {
  grList[isNULL]
  resList <- resList[!isNULL]
  resNames <- names(grList)[!isNULL]
  pLabels <- pLabels[!isNULL]
}

resL <- sapply(resList, length)
graphNames <- rep(resNames, times=resL)
pathwayNames <- rep(pLabels, times=resL)

graphList <- NULL
for (res in resList) {
  grl <- lapply(res, FUN=function(x) {
    x$graph
  })
  graphList <- c(graphList, as.list(grl))
}

ndims <- NULL
for (res in resList) {
  ndim <- sapply(res, FUN=function(x) {
    x$k
  })
  ndims <- c(ndims, ndim)
}

pKEGG <- NULL
for (res in resList) {
  pp <- sapply(res, FUN=function(x) {
    x$p.value
  })
  pKEGG <- cbind(pKEGG, pp)
}
colnames(pKEGG) <- graphNames
rn <- rownames(pKEGG)
rownames(pKEGG)[grep("Fourier", rn)] <- paste("T2 (", round(100*prop), "% Fourier components)", sep="")

if (exists("maPalette", mode="function")) {
  pal <- maPalette(low="red", high="green", mid="black", k=100)
} else {
  pal <- heat.colors(100)
}
shift <- tts # Plot t-statistics
names(shift) <- translateGeneID2KEGGID(names(tts))
fSignif <- which(pKEGG[2,] < 0.05)
fSignif <- fSignif[order(pKEGG[2,fSignif])]


###################################################
### code chunk number 6: path1
###################################################
gIdx <- fSignif[1]
gr <- graphList[[gIdx]]
mm <- match(translateKEGGID2GeneID(nodes(gr)), rownames(annData))
dn <- annData[mm, "NCBI.gene.symbol"]
res <- plotValuedGraph(gr, values=shift, nodeLabels=dn, qMax=0.95, colorPalette=pal, height=40, lwd=1, cex=0.3, verbose=verbose)
stext(side=3, pos=0, pathwayNames[gIdx])
ps <- signif(pKEGG[, gIdx],2)
txt1 <- paste("p(T2)=", ps[1], sep="")
txt2 <- paste("p(T2F[", ndims[gIdx], "])=", ps[2], sep="")
txt <- paste(txt1, txt2, sep="\n")
stext(side=3, pos=1, txt)
image.plot(legend.only=TRUE, zlim=range(res$breaks), col=pal, legend.shrink=0.3, legend.width=0.8, legend.lab="t-scores", legend.mar=3.3) 


###################################################
### code chunk number 7: path2
###################################################
gIdx <- fSignif[5]
gr <- graphList[[gIdx]]
mm <- match(translateKEGGID2GeneID(nodes(gr)), rownames(annData))
dn <- annData[mm, "NCBI.gene.symbol"]
res <- plotValuedGraph(gr, values=shift, nodeLabels=dn, qMax=0.95, colorPalette=pal, height=40, lwd=1, cex=0.3, verbose=verbose)
stext(side=3, pos=0, pathwayNames[gIdx])
ps <- signif(pKEGG[, gIdx],2)
txt1 <- paste("p(T2)=", ps[1], sep="")
txt2 <- paste("p(T2F[", ndims[gIdx], "])=", ps[2], sep="")
txt <- paste(txt1, txt2, sep="\n")
stext(side=3, pos=1, txt)
image.plot(legend.only=TRUE, zlim=range(res$breaks), col=pal, legend.shrink=0.3, legend.width=0.8, legend.lab="t-scores", legend.mar=3.3) 


###################################################
### code chunk number 8: showpath1
###################################################
gIdx <- fSignif[1]
gr <- graphList[[gIdx]]
mm <- match(translateKEGGID2GeneID(nodes(gr)), rownames(annData))
dn <- annData[mm, "NCBI.gene.symbol"]
res <- plotValuedGraph(gr, values=shift, nodeLabels=dn, qMax=0.95, colorPalette=pal, height=40, lwd=1, cex=0.3, verbose=verbose)
stext(side=3, pos=0, pathwayNames[gIdx])
ps <- signif(pKEGG[, gIdx],2)
txt1 <- paste("p(T2)=", ps[1], sep="")
txt2 <- paste("p(T2F[", ndims[gIdx], "])=", ps[2], sep="")
txt <- paste(txt1, txt2, sep="\n")
stext(side=3, pos=1, txt)
image.plot(legend.only=TRUE, zlim=range(res$breaks), col=pal, legend.shrink=0.3, legend.width=0.8, legend.lab="t-scores", legend.mar=3.3) 


###################################################
### code chunk number 9: showpath2
###################################################
gIdx <- fSignif[5]
gr <- graphList[[gIdx]]
mm <- match(translateKEGGID2GeneID(nodes(gr)), rownames(annData))
dn <- annData[mm, "NCBI.gene.symbol"]
res <- plotValuedGraph(gr, values=shift, nodeLabels=dn, qMax=0.95, colorPalette=pal, height=40, lwd=1, cex=0.3, verbose=verbose)
stext(side=3, pos=0, pathwayNames[gIdx])
ps <- signif(pKEGG[, gIdx],2)
txt1 <- paste("p(T2)=", ps[1], sep="")
txt2 <- paste("p(T2F[", ndims[gIdx], "])=", ps[2], sep="")
txt <- paste(txt1, txt2, sep="\n")
stext(side=3, pos=1, txt)
image.plot(legend.only=TRUE, zlim=range(res$breaks), col=pal, legend.shrink=0.3, legend.width=0.8, legend.lab="t-scores", legend.mar=3.3)