### R code from vignette source 'vignettes/predictionet/inst/doc/predictionet.Rnw'

###################################################
### code chunk number 1: run
###################################################
% Benjamin Haibe-Kains
%
% 
% October 5, 2011

\documentclass[a4paper,11pt]{article}

%\VignetteIndexEntry{predictionet}
%\VignetteDepends{predictionet}
%\VignetteKeywords{predictionet}
%\VignettePackage{predictionet}

\usepackage{helvet}
\usepackage[helvet]{sfmath}
\renewcommand{\familydefault}{\sfdefault}
\usepackage[american]{babel}
\usepackage{url}
\usepackage{hyperref}
\usepackage{a4wide}
%\usepackage{colortbl}
%\usepackage{longtable}
%\usepackage{multirow}
%\usepackage{rotating}
\usepackage{natbib}
\usepackage{graphicx}
\usepackage{authblk}
%additional packages
%\usepackage{minitoc}
%\usepackage{amsmath}
%\usepackage{bm}
%\usepackage{algorithm}
%\usepackage{algpseudocode}
\usepackage{wasysym}
\usepackage{tikz}

\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}

\title{\Rpackage{predictionet}: a package for inferring predictive networks from high-dimensional genomic data}
\author[1,2]{Benjamin Haibe-Kains}
\author[3]{Catharina Olsen}
\author[3]{Gianluca Bontempi}
\author[1,2]{John Quackenbush}
\affil[1]{Computational Biology and Functional Genomics Laboratory, Dana-Farber Cancer Institute, Harvard School of Public Health}
\affil[2]{Center for Cancer Computational Biology, Dana-Farber Cancer Institute}
\affil[3]{Machine Learning Group, Universit\'{e} Libre de Bruxelles}

\begin{document}

\SweaveOpts{highlight=TRUE, tidy=TRUE, keep.space=TRUE, keep.blank.space=FALSE, keep.comment=TRUE, keep.source=TRUE}



###################################################
### code chunk number 2: setup
###################################################
if(dorun) {
## initialize seed
set.seed(123)
}


###################################################
### code chunk number 3: loadlib
###################################################
if(dorun) {
library(predictionet)
}


###################################################
### code chunk number 4: preparepriors
###################################################
if(dorun) {
## RAS-related genes
genes.ras <- colnames(data.ras)

## read priors generated by the Predictive Networks web application
pn.priors <- read.csv(system.file("extdata/priors_ras_bild2006_pnwebapp.csv", package="predictionet"), stringsAsFactors=FALSE)
## the column names should be: subject, predicate, object, direction, evidence, sentence, article, network

## remove special characters in the gene symbols
pn.priors[ ,"subject"] <- gsub(pattern="[-]|[+]|[*]|[%]|[$]|[#]|[{]|[}]|[[]|[]]|[|]|[\\^]", replacement=".", x=pn.priors[ ,"subject"])
pn.priors[ ,"object"] <- gsub(pattern="[-]|[+]|[*]|[%]|[$]|[#]|[{]|[}]|[[]|[]]|[|]|[\\^]", replacement=".", x=pn.priors[ ,"object"])
genes.ras <- gsub(pattern="[-]|[+]|[*]|[%]|[$]|[#]|[{]|[}]|[[]|[]]|[|]|[\\^]", replacement=".", x=genes.ras)

## missing values
pn.priors[!is.na(pn.priors) & (pn.priors == "" | pn.priors == " " | pn.priors == "N/A")] <- NA

## select only the interactions in which the genes are comprised in our gene expression dataset
myx <- is.element(pn.priors[ , "subject"], genes.ras) & is.element(pn.priors[ , "object"], genes.ras)
pn.priors <- pn.priors[myx, , drop=FALSE]
}


###################################################
### code chunk number 5: priorsrandom
###################################################
if(dorun) {
pn.priors <- pn.priors[sample(1:nrow(pn.priors)), ]
}


###################################################
### code chunk number 6: tablepriors
###################################################
if(dorun) {
## print the first rows of 'pn.priors'
print(head(pn.priors))
}


###################################################
### code chunk number 7: preparepriorscount
###################################################
if(dorun) {
## build prior counts
pn.priors.counts <- matrix(0, nrow=length(genes.ras), ncol=length(genes.ras), dimnames=list(genes.ras, genes.ras))

for(i in 1:nrow(pn.priors)) {
	switch(tolower(pn.priors[i, "direction"]), 
	"right"={ pn.priors.counts[pn.priors[i, "subject"], pn.priors[i, "object"]] <- pn.priors.counts[pn.priors[i, "subject"], pn.priors[i, "object"]] + ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0) }, 
	"left"={ pn.priors.counts[pn.priors[i, "object"], pn.priors[i, "subject"]] <- pn.priors.counts[pn.priors[i, "object"], pn.priors[i, "subject"]] +  ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0) }, 
	{ pn.priors.counts[pn.priors[i, "subject"], pn.priors[i, "object"]] <- pn.priors.counts[pn.priors[i, "subject"], pn.priors[i, "object"]] +  ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0)
	if(pn.priors[i, "object"] != pn.priors[i, "subject"]) { pn.priors.counts[pn.priors[i, "object"], pn.priors[i, "subject"]] <- pn.priors.counts[pn.priors[i, "object"], pn.priors[i, "subject"]] +  ifelse(!is.na(pn.priors[i, "evidence"]), ifelse(tolower(pn.priors[i, "evidence"]) == "positive", +1, -1), 0) } })
}
## negative count represent evidence for ABSENCE of an interaction, positive otherwise
}


###################################################
### code chunk number 8: priorscounttab
###################################################
if(dorun) {
print(table(pn.priors.counts))
}


###################################################
### code chunk number 9: foopredictionet (eval = FALSE)
###################################################
## if(dorun) {
## library(help=predictionet)
## }


###################################################
### code chunk number 10: loadpredictionet (eval = FALSE)
###################################################
## if(dorun) {
## help(expO.colon.ras)
## help(jorissen.colon.ras)
## }


###################################################
### code chunk number 11: helpnetinf (eval = FALSE)
###################################################
## if(dorun) {
## help(netinf)
## }


###################################################
### code chunk number 12: firstnetinf
###################################################
if(dorun) {
## number of genes to select for the analysis
genen <- 30
## select only the top genes
goi <- dimnames(annot.ras)[[1]][order(abs(log2(annot.ras[ ,"fold.change"])), decreasing=TRUE)[1:genen]]
}


###################################################
### code chunk number 13: firstnetinf2
###################################################
if(dorun) {
mynet <- netinf(data=data.ras[ ,goi], priors=pn.priors.counts[goi,goi], priors.count=TRUE, priors.weight=0.5, maxparents=4, method="regrnet", seed=54321)
}


###################################################
### code chunk number 14: regrnetopofig
###################################################
if(dorun) {
## network topology
#mytopoglobal <- net2topo(net=mynet)
mytopoglobal <- mynet$topology
library(network)
xnet <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
mynetlayout <- plot.network(x=xnet, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.6)
}


###################################################
### code chunk number 15: edgecoldiff
###################################################
if(dorun) {
## preparing colors
mycols <- c("blue", "green", "red")
names(mycols) <- c("prior", "data", "both")
myedgecols <- matrix("white", nrow=nrow(mytopoglobal), ncol=ncol(mytopoglobal), dimnames=dimnames(mytopoglobal))
## prior only
myedgecols[mytopoglobal == 0 & pn.priors.counts[goi,goi] >= 1] <- mycols["prior"]
## data only
myedgecols[mytopoglobal == 1 & pn.priors.counts[goi,goi] < 1] <- mycols["data"]
## both in priors and data
myedgecols[mytopoglobal == 1 & pn.priors.counts[goi,goi] >= 1] <- mycols["both"]
}


###################################################
### code chunk number 16: edgecoldiffig
###################################################
if(dorun) {
mytopoglobal2 <- (myedgecols != "white") + 0
## network topology
xnet2 <- network(x=mytopoglobal2, matrix.type="adjacency", directed=TRUE, loops=TRUE, vertex.attrnames=dimnames(mytopoglobal2)[[1]])
plot.network(x=xnet2, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=myedgecols, coord=mynetlayout)
}


###################################################
### code chunk number 17: regrnetinfcv
###################################################
if(dorun) {
myres.cv <- netinf.cv(data=data.ras[ ,goi], categories=3, priors=pn.priors.counts[goi,goi], maxparents=4, priors.weight=0.5, method="regrnet", nfold=10, seed=54321)
}


###################################################
### code chunk number 18: rescv
###################################################
if(dorun) {
print(str(myres.cv, 1))
}


###################################################
### code chunk number 19: edgecol
###################################################
if(dorun) {
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- as.character(ii/10)
myedgecols <- matrix("#00000000", nrow=nrow(mytopoglobal), ncol=ncol(mytopoglobal), dimnames=dimnames(mytopoglobal))
for(k in 1:length(mycols)) { myedgecols[myres.cv$edge.stability == names(mycols)[k]] <- mycols[k] }
myedgecols[!mytopoglobal] <- "#00000000"
}


###################################################
### code chunk number 20: edgestabfig0 (eval = FALSE)
###################################################
## if(dorun) {
## def.par <- par(no.readonly=TRUE)
## layout(rbind(1,2), heights=rbind(8,1))
## par(mar=c(1,1,1,1))
## ## network topology
## xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=myedgecols, coord=mynetlayout)
## par(mar=c(0,3,1,3))
## plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="Stability scale", cex.main=1)
## rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
## text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
## par(def.par)
## }


###################################################
### code chunk number 21: edgestabfig
###################################################
if(dorun) {
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=myedgecols, coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="Stability scale", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)
}


###################################################
### code chunk number 22: genecolr2
###################################################
if(dorun) {
myr2 <- apply(myres.cv$prediction.score.cv$r2, 2, mean, na.rm=TRUE)
myr2 <- round(myr2*10)/10
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- as.character(ii/10)
myvertexcols <- mycols[match(myr2, names(mycols))]
names(myvertexcols) <- names(myr2)
}


###################################################
### code chunk number 23: genepredabr2fig0 (eval = FALSE)
###################################################
## if(dorun) {
## def.par <- par(no.readonly=TRUE)
## layout(rbind(1,2), heights=rbind(8,1))
## par(mar=c(1,1,1,1))
## ## network topology
## xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
## par(mar=c(0,3,1,3))
## plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$R^2$", cex.main=1)
## rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
## text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
## par(def.par)
## }


###################################################
### code chunk number 24: genepredabr2fig
###################################################
if(dorun) {
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$R^2$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)
}


###################################################
### code chunk number 25: genecolmcc
###################################################
if(dorun) {
mymcc <- apply(myres.cv$prediction.score.cv$mcc, 2, mean, na.rm=TRUE)
mymcc <- round(mymcc*20)/20
mymcc[mymcc < 0.5] <- 0.5
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- seq(0.5, 1, 0.05)
myvertexcols <- mycols[match(mymcc, names(mycols))]
names(myvertexcols) <- names(mymcc)
}


###################################################
### code chunk number 26: genepredabmccfig0 (eval = FALSE)
###################################################
## if(dorun) {
## def.par <- par(no.readonly=TRUE)
## layout(rbind(1,2), heights=rbind(8,1))
## par(mar=c(1,1,1,1))
## ## network topology
## xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
## par(mar=c(0,3,1,3))
## plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$MCC$", cex.main=1)
## rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
## text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
## par(def.par)
## }


###################################################
### code chunk number 27: genepredabmccfig
###################################################
if(dorun) {
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$MCC$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)
}


###################################################
### code chunk number 28: validationpred
###################################################
if(dorun) {
## compute predictions
mynet.test.pred <-  netinf.predict(net=mynet, data=data2.ras[ ,goi])
## performance estimation: R2
mynet.test.r2 <- pred.score(data=data2.ras[ ,goi], pred=mynet.test.pred, method="r2")
## performance estimation: MCC
mynet.test.mcc <- pred.score(data=data2.ras[ ,goi], categories=3, pred=mynet.test.pred, method="mcc")
}


###################################################
### code chunk number 29: genecolr2test
###################################################
if(dorun) {
mynet.test.r2 <- round(mynet.test.r2*10)/10
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- as.character(ii/10)
myvertexcols <- mycols[match(mynet.test.r2, names(mycols))]
names(myvertexcols) <- names(mynet.test.r2)
}


###################################################
### code chunk number 30: genepredabr2testfig0 (eval = FALSE)
###################################################
## if(dorun) {
## def.par <- par(no.readonly=TRUE)
## layout(rbind(1,2), heights=rbind(8,1))
## par(mar=c(1,1,1,1))
## ## network topology
## xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
## par(mar=c(0,3,1,3))
## plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$R^2$", cex.main=1)
## rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
## text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
## par(def.par)
## }


###################################################
### code chunk number 31: genepredabr2testfig
###################################################
if(dorun) {
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$R^2$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)
}


###################################################
### code chunk number 32: genecolmcctest
###################################################
if(dorun) {
mynet.test.mcc <- round(mynet.test.mcc*20)/20
mynet.test.mcc[mynet.test.mcc < 0.5] <- 0.5
## preparing colors
ii <- 0:10
mycols <- c("#BEBEBE", rev(rainbow(10, v=0.8, alpha=0.5)))
names(mycols) <- seq(0.5, 1, 0.05)
myvertexcols <- mycols[match(mynet.test.mcc, names(mycols))]
names(myvertexcols) <- names(mynet.test.mcc)
}


###################################################
### code chunk number 33: genepredabmcctestfig0 (eval = FALSE)
###################################################
## if(dorun) {
## def.par <- par(no.readonly=TRUE)
## layout(rbind(1,2), heights=rbind(8,1))
## par(mar=c(1,1,1,1))
## ## network topology
## xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
## plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
## par(mar=c(0,3,1,3))
## plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$MCC$", cex.main=1)
## rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
## text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
## par(def.par)
## }


###################################################
### code chunk number 34: genepredabmcctestfig
###################################################
if(dorun) {
def.par <- par(no.readonly=TRUE)
layout(rbind(1,2), heights=rbind(8,1))
par(mar=c(1,1,1,1))
## network topology
xnet3 <- network(x=mytopoglobal, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopoglobal)[[1]])
plot.network(x=xnet3, displayisolates=TRUE, displaylabels=TRUE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col=myvertexcols, jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout)
par(mar=c(0,3,1,3))
plot(ii+1, ii+10/6+1, bty="n", type="n", yaxt="n", xaxt="n", ylab="", xlab ="", main="$MCC$", cex.main=1)
rect(xleft=ii+0.5, ybottom=3, xright=ii+1.4, ytop=10+3, col=mycols, border="grey")
text(ii+1, y=2.4, labels=names(mycols), pos=3, cex=1)
par(def.par)
}


###################################################
### code chunk number 35: regnetcvexport (eval = FALSE)
###################################################
## if(dorun) {
## netinf2gml(object=myres.cv, file="regrnet_expO_cv")
## }


###################################################
### code chunk number 36: regnetcvpriorsweight
###################################################
if(dorun) {
## priors.weight 0
myres.cv.pw0 <- netinf.cv(data=data.ras[ ,goi], categories=3, priors=pn.priors.counts[goi,goi], maxparents=4, priors.weight=0, method="regrnet", nfold=10, seed=54321)
## priors.weight 0.25
myres.cv.pw025 <- netinf.cv(data=data.ras[ ,goi], categories=3, priors=pn.priors.counts[goi,goi], maxparents=4, priors.weight=0.25, method="regrnet", nfold=10, seed=54321)
## priors.weight 0.5
myres.cv.pw050 <- myres.cv
## priors.weight 0.75
myres.cv.pw075 <- netinf.cv(data=data.ras[ ,goi], categories=3, priors=pn.priors.counts[goi,goi], maxparents=4, priors.weight=0.75, method="regrnet", nfold=10, seed=54321)
## priors.weight 0
myres.cv.pw1 <- netinf.cv(data=data.ras[ ,goi], categories=3, priors=pn.priors.counts[goi,goi], maxparents=4, priors.weight=1, method="regrnet", nfold=10, seed=54321)
}


###################################################
### code chunk number 37: regnetcvpriorsweightfig20 (eval = FALSE)
###################################################
## if(dorun) {
## def.par <- par(no.readonly=TRUE)
## layout(mat=matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
## ## priors.weight 0
## mytopot <- myres.cv.pw0$topology
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0 (data only)")
## ## priors.weight 0.25
## mytopot <- myres.cv.pw025$topology
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.25")
## ## priors.weight 0.75
## mytopot <- myres.cv.pw075$topology
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.75")
## ## priors.weight 1
## mytopot <-  myres.cv.pw1$topology
## xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 1 (priors only)")
## par(def.par)
## }


###################################################
### code chunk number 38: regnetcvpriorsweightfig2
###################################################
if(dorun) {
def.par <- par(no.readonly=TRUE)
layout(mat=matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
## priors.weight 0
mytopot <- myres.cv.pw0$topology
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0 (data only)")
## priors.weight 0.25
mytopot <- myres.cv.pw025$topology
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.25")
## priors.weight 0.75
mytopot <- myres.cv.pw075$topology
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 0.75")
## priors.weight 1
mytopot <-  myres.cv.pw1$topology
xnett <- network(x=mytopot, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col="black", coord=mynetlayout, main="Priors weight = 1 (priors only)")
par(def.par)
}


###################################################
### code chunk number 39: boxplotstabpwcvfig2
###################################################
if(dorun) {
gg <- c("0", "0.25", "0.50", "0.75", "1")
mystab.cv.pw <- list(myres.cv.pw0$edge.stability[myres.cv.pw0$topology == 1], myres.cv.pw025$edge.stability[myres.cv.pw025$topology == 1], myres.cv.pw050$edge.stability[myres.cv.pw050$topology == 1], myres.cv.pw075$edge.stability[myres.cv.pw075$topology == 1], myres.cv.pw1$edge.stability[myres.cv.pw1$topology == 1])
names(mystab.cv.pw) <- gg
boxplot(mystab.cv.pw, xlab="priors.weight", ylim=c(0, 1), ylab="Edge stability", border="grey", col="lightgrey", outline=FALSE)
points(x=jitter(x=rep(1:length(mystab.cv.pw), times=sapply(mystab.cv.pw, length)), amount=0.15), y=unlist(mystab.cv.pw), cex=0.55, pch=16, col="royalblue")
}


###################################################
### code chunk number 40: boxplotr2pwcvfig2
###################################################
if(dorun) {
gg <- c("0", "0.25", "0.50", "0.75", "1")
myr2.cv.pw <- cbind(apply(myres.cv.pw0$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw025$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw050$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw075$prediction.score.cv$r2, 2, mean, na.rm=TRUE), apply(myres.cv.pw1$prediction.score.cv$r2, 2, mean, na.rm=TRUE))
colnames(myr2.cv.pw) <- gg
gg <- factor(rep(gg, each=genen), levels=gg)
boxplot(myr2.cv.pw, xlab="priors.weight", ylim=c(0, 1), ylab="$R^2$", border="grey", col="lightgrey", outline=FALSE)
points(x=jitter(x=rep(1:ncol(myr2.cv.pw), times=nrow(myr2.cv.pw)), amount=0.15), y=as.vector(t(myr2.cv.pw)), cex=0.55, pch=16, col="royalblue")
}


###################################################
### code chunk number 41: regnetcvpriorsweightcompr2
###################################################
if(dorun) {
## Friedman test to test whether at least one of the networks gives statstically different predictive ability
print(friedman.test(y=myr2.cv.pw))

## Pairwise Wilcoxon Rank Sum tests
print(pairwise.wilcox.test(x=as.vector(myr2.cv.pw), g=gg, paired=TRUE, exact=FALSE, alternative="two.sided", p.adjust.method="holm"))
}


###################################################
### code chunk number 42: netinfcv12
###################################################
if(dorun) {
## expO
myres21.cv <- netinf.cv(data=data.ras[ ,goi], categories=3, priors=priors2.ras[goi,goi], maxparents=4, priors.weight=0, method="regrnet", nfold=10, seed=54321)

## jorissen
myres22.cv <- netinf.cv(data=data2.ras[ ,goi], categories=3, priors=priors2.ras[goi,goi], maxparents=4, priors.weight=0, method="regrnet", nfold=10, seed=54321)
}


###################################################
### code chunk number 43: regnetcvtopo1topo2fig20 (eval = FALSE)
###################################################
## if(dorun) {
## topo1 <- myres21.cv$topology
## topo2 <- myres22.cv$topology
## ## preparing colors
## myedgecols <- matrix("white", nrow=nrow(topo1), ncol=ncol(topo1), dimnames=dimnames(topo1))
## myedgecols[topo1== 1 & topo2 == 1 & pn.priors.counts[rownames(topo1), colnames(topo1)] > 0] <- "orange"
## myedgecols[topo1== 1 & topo2 == 1 & pn.priors.counts[rownames(topo1), colnames(topo1)]  <= 0] <- "red"
## 
## def.par <- par(no.readonly=TRUE)
## layout(mat=matrix(1:2, nrow=1, ncol=2, byrow=TRUE))
## ## expO
## mycolt <- myedgecols
## mycolt[myedgecols == "white" & topo1 == 1 ] <- "gray"
## xnett <- network(x=topo1, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=mycolt, coord=mynetlayout, main="Dataset: expO.colon.ras")
## ## jorissen
## mycolt <- myedgecols
## mycolt[myedgecols == "white" & topo2 == 1 ] <- "gray"
## xnett <- network(x=topo2, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
## plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=mycolt, coord=mynetlayout, main="Dataset: jorissen.colon.ras")
## par(def.par)
## }


###################################################
### code chunk number 44: regnetcvtopo1topo2fig2
###################################################
if(dorun) {
topo1 <- myres21.cv$topology
topo2 <- myres22.cv$topology
## preparing colors
myedgecols <- matrix("white", nrow=nrow(topo1), ncol=ncol(topo1), dimnames=dimnames(topo1))
myedgecols[topo1== 1 & topo2 == 1 & pn.priors.counts[rownames(topo1), colnames(topo1)] > 0] <- "orange"
myedgecols[topo1== 1 & topo2 == 1 & pn.priors.counts[rownames(topo1), colnames(topo1)]  <= 0] <- "red"

def.par <- par(no.readonly=TRUE)
layout(mat=matrix(1:2, nrow=1, ncol=2, byrow=TRUE))
## expO
mycolt <- myedgecols
mycolt[myedgecols == "white" & topo1 == 1 ] <- "gray"
xnett <- network(x=topo1, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=mycolt, coord=mynetlayout, main="Dataset: expO.colon.ras")
## jorissen
mycolt <- myedgecols
mycolt[myedgecols == "white" & topo2 == 1 ] <- "gray"
xnett <- network(x=topo2, matrix.type="adjacency", directed=TRUE, loops=FALSE, vertex.attrnames=dimnames(mytopot)[[1]])
plot.network(x=xnett, displayisolates=TRUE, displaylabels=FALSE, boxed.labels=FALSE, label.pos=0, arrowhead.cex=1.5, vertex.cex=2, vertex.col="royalblue", jitter=FALSE, pad=0.5, edge.col=mycolt, coord=mynetlayout, main="Dataset: jorissen.colon.ras")
par(def.par)
}


###################################################
### code chunk number 45: regnetcvtopo1topo2stabfig20 (eval = FALSE)
###################################################
## if(dorun) {
## def.par <- par(no.readonly=TRUE)
## layout(mat=matrix(1:2, nrow=1, ncol=2, byrow=TRUE))
## ## expO
## stab2 <- list("specific"=myres21.cv$edge.stability[(topo1 == 1 & topo2 == 0)], "common"=myres21.cv$edge.stability[topo1 == 1 & topo2 == 1])
## wt <- wilcox.test(x=stab2$specific, y=stab2$common)
## boxplot(stab2, ylab="Stability", ylim=c(0, 1), xlab="", border="grey", col="lightgrey", outline=FALSE, sub=sprintf("Wilcoxon test p-value = %.1E", wt$p.value), main="Dataset: expO.colon.ras")
## points(x=jitter(x=rep(1:length(stab2), times=sapply(stab2, length)), amount=0.15), y=unlist(stab2), cex=0.55, pch=16, col="royalblue")
## ## jorissen
## stab2 <- list("specific"=myres22.cv$edge.stability[(topo1 == 0 & topo2 == 1)], "common"=myres22.cv$edge.stability[topo1 == 1 & topo2 == 1])
## wt <- wilcox.test(x=stab2$specific, y=stab2$common)
## boxplot(stab2, ylab="Stability", ylim=c(0, 1), xlab="", border="grey", col="lightgrey", outline=FALSE, sub=sprintf("Wilcoxon test p-value = %.1E", wt$p.value), main="Dataset: jorissen.colon.ras")
## points(x=jitter(x=rep(1:length(stab2), times=sapply(stab2, length)), amount=0.15), y=unlist(stab2), cex=0.55, pch=16, col="royalblue")
## par(def.par)
## }


###################################################
### code chunk number 46: regnetcvtopo1topo2stabfig2
###################################################
if(dorun) {
def.par <- par(no.readonly=TRUE)
layout(mat=matrix(1:2, nrow=1, ncol=2, byrow=TRUE))
## expO
stab2 <- list("specific"=myres21.cv$edge.stability[(topo1 == 1 & topo2 == 0)], "common"=myres21.cv$edge.stability[topo1 == 1 & topo2 == 1])
wt <- wilcox.test(x=stab2$specific, y=stab2$common)
boxplot(stab2, ylab="Stability", ylim=c(0, 1), xlab="", border="grey", col="lightgrey", outline=FALSE, sub=sprintf("Wilcoxon test p-value = %.1E", wt$p.value), main="Dataset: expO.colon.ras")
points(x=jitter(x=rep(1:length(stab2), times=sapply(stab2, length)), amount=0.15), y=unlist(stab2), cex=0.55, pch=16, col="royalblue")
## jorissen
stab2 <- list("specific"=myres22.cv$edge.stability[(topo1 == 0 & topo2 == 1)], "common"=myres22.cv$edge.stability[topo1 == 1 & topo2 == 1])
wt <- wilcox.test(x=stab2$specific, y=stab2$common)
boxplot(stab2, ylab="Stability", ylim=c(0, 1), xlab="", border="grey", col="lightgrey", outline=FALSE, sub=sprintf("Wilcoxon test p-value = %.1E", wt$p.value), main="Dataset: jorissen.colon.ras")
points(x=jitter(x=rep(1:length(stab2), times=sapply(stab2, length)), amount=0.15), y=unlist(stab2), cex=0.55, pch=16, col="royalblue")
par(def.par)
}


###################################################
### code chunk number 47: regnetcvtopo1topo2predfig2
###################################################
if(dorun) {
pred2 <- list("expO"=apply(myres21.cv$prediction.score.cv$r2, 2, mean, na.rm=TRUE), "jorissen"=apply(myres22.cv$prediction.score.cv$r2, 2, mean, na.rm=TRUE))
plot(x=pred2$expO, y=pred2$jorissen, xlim=c(0, 0.6), ylim=c(0, 0.6), pch=16, col="royalblue", xlab="Prediction scores in expO.colon.ras", ylab="Prediction scores in jorissen.colon.ras")
legend(x="bottomright", legend=sprintf("cor = %.2g", cor(pred2$expO, pred2$jorissen, method="spearman", use="complete.obs")), bty="n")
}


###################################################
### code chunk number 48: sessioninfo
###################################################
if(dorun) {
toLatex(sessionInfo(), locale=FALSE)
}