@@ -1,7 +1,7 @@

 Package: DEGraph

 Title: Two-sample tests on a graph

-Version: 1.1.0

-Date: 2010-10-14

+Version: 1.2.0

+Date: 2011-03-23

 Author: Laurent Jacob, Pierre Neuvial and Sandrine Dudoit

 Maintainer: Laurent Jacob <laurent.jacob@gmail.com>

 Description: DEGraph implements recent hypothesis testing methods

@@ -16,7 +16,7 @@ Description: DEGraph implements recent hypothesis testing methods

   and tools to visualize the results.

 License: GPL-3

 LazyLoad: yes

-Imports: graph, KEGGgraph, lattice, mvtnorm, R.methodsS3, RBGL, Rgraphviz, rrcov

-Suggests: corpcor, fields, graph, KEGGgraph, lattice, marray, RBGL, rrcov, Rgraphviz

+Imports: graph, KEGGgraph, lattice, mvtnorm, NCIgraph, R.methodsS3, RBGL, Rgraphviz, rrcov

+Suggests: corpcor, fields, graph, KEGGgraph, lattice, marray, NCIgraph, RBGL, rrcov, Rgraphviz

 Depends: R (>= 2.10.0), R.utils

 biocViews: Microarray, Bioinformatics, DifferentialExpression, GraphsAndNetworks

@@ -25,4 +25,4 @@ importFrom(rrcov, T2.test)

 importFrom(mvtnorm, rmvnorm)

 importFrom(RBGL,connectedComp)

 importFrom(lattice, level.colors)

-

+importFrom(NCIgraph, is.NCIgraph)

@@ -88,11 +88,14 @@ getSignedGraph <- function(graph, positiveInteractionLabels=c("activation", "exp

     cat <- R.utils::cat

     pushState(verbose)

     on.exit(popState(verbose))

-  } 

+  }

-  ## retrieve edge types

-  st <- getSubtype(graph)  ## FIXME: does not work for undirected graphs

+  ## retrieve edge types (method detects whether or not graph is a NCIgraph)

+  ##  if(is.NCIgraph(graph))

+  ##    st <- getSubtype.NCIgraph(graph)

+  ##  else

+  st <- getSubtype(graph) ## FIXME: does not work for undirected graphs

   edgeNames <- names(st)

   intNames <- sapply(st, FUN=function(x) {

     x$subtype@name

@@ -134,7 +137,7 @@ getSignedGraph <- function(graph, positiveInteractionLabels=c("activation", "exp

     etrNames <- edgeNames[-c(idxPos, idxNeg)]

     ftNames <- sapply(etrNames, FUN=function(edge) {

-      unlist(strsplit(edge, "~"))

+      unlist(strsplit(edge, "~|\\|"))

     })

     oldGraph <- graph

@@ -161,7 +164,7 @@ getSignedGraph <- function(graph, positiveInteractionLabels=c("activation", "exp

   if (length(negNames)) {

     ## set negative edges to -1

     ftNames <- sapply(negNames, FUN=function(edge) {

-      unlist(strsplit(edge, "~"))

+      unlist(strsplit(edge, "~|\\|"))

     })

     nodeNames <- colnames(signMat)

     for (ii in seq(length=ncol(ftNames))) {

@@ -180,7 +183,7 @@ getSignedGraph <- function(graph, positiveInteractionLabels=c("activation", "exp

   if (length(posNames)) {

     ## check that positive edges are 1

     ftNames <- sapply(posNames, FUN=function(edge) {

-      unlist(strsplit(edge, "~"))

+      unlist(strsplit(edge, "~|\\|"))

     })

     nodeNames <- colnames(signMat)

     for (ii in seq(length=ncol(ftNames))) {

@@ -201,6 +204,8 @@ getSignedGraph <- function(graph, positiveInteractionLabels=c("activation", "exp

 ############################################################################

 ## HISTORY:

+## 2011-03-06

+## o Dealing with NCI graphs

 ## 2010-10-08

 ## o Now validating argument 'verbose'.

 ## o Updated arguments.

@@ -71,6 +71,7 @@ graph.T2.test <- function(X1, X2, G=NULL, lfA=NULL, ..., k=ncol(X1))

     }    

   U <- lfA$U

+  ## print(U)

   egVal <- lfA$l

   kIdx <- (egVal <= max(egVal[k],tol)) #lfA$kIdx

   rk <- max(which(kIdx)) ## "round" the number of kept eigenvectors to take into account eigenvalue multiplicity

@@ -179,11 +179,26 @@ plotValuedGraph <- function(graph, values=NULL, nodeLabels=nodes(graph), qMax=0.

     if (symmetrizeArrows) {

       edgeRenderInfo(graph) <- list(arrowtail=eArrowhead)

     }

-  } else {

-    ## graphAM

-    if (inherits(graph, "graphAM")) {

-     ahd <- rep("normal", length(ed))

-     names(ahd) <- names(ed)

+  }

+  else

+    if(is.NCIgraph(graph)) ## NCIgraph

+      {

+        eTypeDictionnary <- c('normal','tee')

+        eColDictionnary <- c('red','blue')

+        names(eTypeDictionnary) <- names(eColDictionnary) <- c('activation','inhibition')

+        eTypes <- unlist(lapply(graph@edgeData@data,FUN=function(e) tolower(e$edgeType)))

+        ERIarrowhead =eTypeDictionnary[eTypes]

+        ERIcol =eColDictionnary[eTypes]

+        eNames <- sapply(names(graph@edgeData@data),FUN=function(e) gsub(e,pattern='\\|',replacement='~'))

+        names(ERIarrowhead) <- names(ERIcol) <- eNames

+        edgeRenderInfo(graph) <- list(arrowhead=ERIarrowhead,col=ERIcol)

+      }

+    else

+      {

+        ## graphAM

+        if (inherits(graph, "graphAM")) {

+          ahd <- rep("normal", length(ed))

+          names(ahd) <- names(ed)

      adjMat <- graph@adjMat

      ## not useful: adjacency matrix is not signed...

    }

@@ -100,7 +100,13 @@ testOneGraph <- function(graph, data, classes, useInteractionSigns=TRUE, ..., ve

   ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   verbose && enter(verbose, "Keeping genes in the graph *and* the expression data set")

   dataGN <- rownames(data)

-  graphGN <- translateKEGGID2GeneID(nodes(graph))

+  if(is.NCIgraph(graph))

+    {

+      graphGN <- translateNCI2GeneID(graph)

+      names(graphGN) <- NULL

+    }

+  else

+    graphGN <- translateKEGGID2GeneID(nodes(graph))

   commonGN <- intersect(dataGN, graphGN)

   mm <- match(graphGN, commonGN)

@@ -129,7 +135,13 @@ testOneGraph <- function(graph, data, classes, useInteractionSigns=TRUE, ..., ve

   data <- data[mm, ]

   ## sanity check (ordering should be the same now)

-  graphGN <- translateKEGGID2GeneID(nodes(graph))

+  if(is.NCIgraph(graph))

+    {

+      graphGN <- translateNCI2GeneID(graph)

+      names(graphGN) <- NULL

+    }

+  else

+    graphGN <- translateKEGGID2GeneID(nodes(graph))

   dataGN <- rownames(data)

   stopifnot(all.equal(dataGN, graphGN))

@@ -180,6 +192,8 @@ testOneGraph <- function(graph, data, classes, useInteractionSigns=TRUE, ..., ve

 ############################################################################

 ## HISTORY:

+## 2011-03-06

+## o Dealing with NCI graphs

 ## 2010-10-08

 ## o Now validating argument 'verbose'.

 ############################################################################

@@ -2,7 +2,8 @@

 ## Test all (locally stored) pathways from KEGG on Loi's data

 ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-library("R.utils")

+library('R.utils')

+library('xtable')

 verbose <- Arguments$getVerbose(-8, timestamp=TRUE)

@@ -14,12 +15,22 @@ sourceDirectory(path)

 path <- system.file("demoScripts", package="DEGraph")

 sourceDirectory(path)

+## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+## get all NCI pathways

+## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+

+library('NCIgraph')

+path <- system.file("downloadScripts", package="NCIgraph")

+sourceDirectory(path)

+load(file.path('rawNCINetworks','NCI-cyList.RData'))

+

+grList <- getNCIPathways(cyList=NCI.cyList, verbose=verbose)$pList

 ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ## get all KEGG pathways

 ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-##grList <- getKEGGPathways(organism="hsa", metaTag="non-metabolic", verbose=verbose)

-grList <- getKEGGPathways(organism="hsa", metaTag="non-metabolic", patt="04060", verbose=verbose)

+## grList <- getKEGGPathways(organism="hsa", metaTag="non-metabolic", verbose=verbose)

+## grList <- getKEGGPathways(organism="hsa", metaTag="non-metabolic", patt="04060", verbose=verbose)

 ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ## test all KEGG pathways

@@ -33,34 +44,51 @@ res <- apply(exprData, 1, FUN=function(x) {

 ttpv <- res["p.value", ]

 tts <- res["statistic", ]

-names(ttpv) <- translateGeneID2KEGGID(names(ttpv))

-names(tts) <- translateGeneID2KEGGID(names(tts))

+if(!is.NCIgraph(grList[[1]]))

+  {

+    names(ttpv) <- translateGeneID2KEGGID(names(ttpv))

+    names(tts) <- translateGeneID2KEGGID(names(tts))

+  }

 prop <- 0.2 ## proportion of the spectrum to be retained for T2-Fourier 

 ## Multivariate tests

 resList <- NULL

-##for (ii in 1:40) {

+## for (ii in 1:5) {

 for (ii in seq(along=grList)) {

   verbose && cat(verbose, ii)

   gr <- grList[[ii]]

-  res <- testOneGraph(gr, exprData, classData, verbose=verbose, prop=prop)

-  if(!is.null(res))

+  if(min(length(nodes(gr)),length(gr@edgeData@data))>0) {

+    res <- testOneGraph(gr, exprData, classData, verbose=verbose, prop=prop)

+  } else {

+    res <- NULL

+  }

+  if(!is.null(res)) {

     res <- lapply(res, FUN=function(x) {

-      if(!is.null(x))

-        {

-          ht <- hyper.test(ttpv, nodes(x$graph), thr=0.001)

+      if(!is.null(x)) {

+          if(is.NCIgraph(x$graph)) {

+            geneIDs <- translateNCI2GeneID(x$graph)

+          } else {

+            geneIDs <- translateKEGGID2GeneID(nodes(x$graph))

+          }

+          ht <- hyper.test(ttpv, geneIDs, thr=0.001)

           x$p.value <- c(x$p.value, ht$p.value)

           names(x$p.value)[length(names(x$p.value))] <- "Hypergeometric test"

           return(x)

+        } else {

+          return(NULL)

         }

-      else

-        return(NULL)

     })

+  }

   resList <- c(resList, list(res))

 }

 resNames <- names(grList)

-pLabels <- sapply(grList, attr, "label")

+

+if(is.NCIgraph(grList[[1]])) {

+  pLabels <- names(grList)

+} else {

+  pLabels <- sapply(grList, attr, "label")

+}

 ## get rid of NULL results (no connected component of size > 1)

 isNULL <- sapply(resList, is.null)

@@ -99,7 +127,8 @@ for (res in resList) {

   pp <- sapply(res, FUN=function(x) {

     x$p.value

   })

-  pKEGG <- cbind(pKEGG, pp)

+  if(length(pp))

+    pKEGG <- cbind(pKEGG, pp)

 }

 colnames(pKEGG) <- graphNames

 rn <- rownames(pKEGG)

@@ -159,7 +188,7 @@ title(sprintf("Number of significant genes\n with the three BY-corrected tests a

 ##pathwayNames[graphNames[order(pHG)[[1]]]]

-oo <- order(pF/pH)[1:5]

+oo <- order(pF/pH)[1:25]

 gIdx <- oo[1]

@@ -173,11 +202,19 @@ if (require(marray)) {

 shift <- tts # Plot t-statistics

 ##shift <- -getMeanShift(exprData, classData) # Plot mean shifts

-dn <- getDisplayName(gr, shortLabel=TRUE)

-mm <- match(translateKEGGID2GeneID(nodes(gr)), rownames(annData))

+## dn <- getDisplayName(gr, shortLabel=TRUE)

+

+if(is.NCIgraph(gr)) {

+  mm <- match(translateNCI2GeneID(gr), rownames(annData))

+  nodes(gr) <- translateNCI2GeneID(gr)

+} else {

+  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, verbose=verbose)

+

+##pdf("Loi2008-leukocyte.pdf",width=10,height=10)

+res <- plotValuedGraph(gr, values=shift, nodeLabels=dn, qMax=0.95, colorPalette=pal, height=40, lwd=1, cex=0.7, verbose=verbose)

+##devDone()

 title(paste(pathwayNames[graphNames[gIdx]],"pF=",as.character(signif(pF[gIdx],4)),"pH=",as.character(signif(pH[gIdx],4))))

 ## mm <- match(translateKEGGID2GeneID(nodes(gr)), names(shift))

@@ -188,18 +225,28 @@ title(paste(pathwayNames[graphNames[gIdx]],"pF=",as.character(signif(pF[gIdx],4)

 ## step.

 fSignif <- which(BHpKEGG[2,]<0.05)

-##fSignif <- fSignif[order(BHpKEGG[1,fSignif],decreasing=TRUE)]

-fSignif <- fSignif[order(BHpKEGG[2,fSignif],decreasing=FALSE)] # Order by increasing pF

+fSignif <- fSignif[order(BHpKEGG[1,fSignif],decreasing=TRUE)]

+##fSignif <- fSignif[order(BHpKEGG[2,fSignif],decreasing=FALSE)] # Order by increasing pF

 ## Plot all graphs

 ##dev.new()

-pdf("MI-smoothGraphs-loi.pdf")

+##printFileName <- "Tech-MI-smoothGraphs-loi-dec"

+printFileName <- "Loi-NCI"

+printStarted <- FALSE

+pdf(paste(printFileName,".pdf",sep=""))

 for (gIdx in fSignif) {

   ## Get graph

   gr <- graphList[[gIdx]]

-  mm <- match(translateKEGGID2GeneID(nodes(gr)), rownames(annData))

+

+  if(is.NCIgraph(gr)) {

+    mm <- match(translateNCI2GeneID(gr), rownames(annData))

+    nodes(gr) <- translateNCI2GeneID(gr)

+  } else {

+    mm <- match(translateKEGGID2GeneID(nodes(gr)), rownames(annData))

+  }

   dn <- annData[mm, "NCBI.gene.symbol"]

+

   ## Try to plot

   tt <- try(res <- plotValuedGraph(gr, values=shift, nodeLabels=dn, qMax=0.95, colorPalette=pal, height=40, lwd=1, cex=0.3, verbose=verbose))

   ## Add legend

@@ -210,7 +257,14 @@ for (gIdx in fSignif) {

     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) 

+    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)

+    ## Now print gene table in a tex file

+    cIdx <- match(nodes(gr),names(shift)) #which(names(shift)%in% nodes(gr))

+    gTable <- cbind(dn,signif(shift[cIdx],2),signif(ttpv[cIdx],2))

+    colnames(gTable) <- c("Gene symbol","t-stat","p-value")

+    xt <- xtable(gTable,caption=sprintf("%s, p(Hotelling)=%g, p(netHotelling)=%g",pathwayNames[gIdx], ps[1], ps[2]))

+    print(xt,file=paste(printFileName,"-Tables.tex",sep=""),tabular.environment='longtable',floating=FALSE, append=printStarted)

+    printStarted <- TRUE

   } else {

     warning("check gIdx=", gIdx)

   }

@@ -221,5 +275,3 @@ devDone()

 ts <- format(Sys.time(), "%Y-%m-%d,%X")

 filename <- sprintf("pKEGG,signed,normalized,%s.rda", ts)

 save(resList, pKEGG, graphList, pathwayNames, ndims, file=filename)

-

-

@@ -1,6 +1,10 @@

 Package: DEGraph

 ===================

+Version: 1.2.0 [2011-03-23]

+o Added support of NCI networks through the NCIgraph package.

+o Updated the man page of the Loi2008 data.

+

 Version: 0.3.3 [2010-10-08]

 o Now validating 'verbose' argument in all functions.

@@ -16,6 +16,20 @@ data("Loi2008_DEGraphVignette")

 dim(exprLoi2008)

 head(exprLoi2008)

 }

+

+\details{

+  The original data set corresponds to data processed by RMA and median-centered

+  as available from the GSE6532 GEO archive:

+  http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6532.

+  

+  These data were summarized from the probe set level to the gene level as 

+  follows. The expression level of a gene was defined as the expression level 

+  of the probe set with largest alignment score among all probe sets mapping 

+  to this gene according to the annotation in GSE6532. When the largest 

+  alignment score was achieved by several probe sets, the median expression 

+  level of those probe sets was taken.

+}

+

 \source{

   Loi et al.,

   \emph{Predicting prognosis using molecular profiling in estrogen