Nothing
inst/doc/predictionet.R
In predictionet: Inference for predictive networks designed for (but not limited to) genomic data
### R code from vignette source 'predictionet.Rnw'
###################################################
### code chunk number 1: setup
###################################################
## initialize seed
set.seed(123)
###################################################
### code chunk number 2: loadlib
###################################################
library(predictionet)
###################################################
### code chunk number 3: preparepriors
###################################################
## RAS-related genes
genes.ras <- colnames(data.ras)
## read priors generated by the Predictive Networks web application
pn.priors <- read.csv(system.file(file.path("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 4: priorsrandom
###################################################
pn.priors <- pn.priors[sample(1:nrow(pn.priors)), ]
###################################################
### code chunk number 5: tablepriors
###################################################
print(head(pn.priors))
###################################################
### code chunk number 6: preparepriorscount
###################################################
## 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 7: priorscounttab
###################################################
print(table(pn.priors.counts))
###################################################
### code chunk number 8: foopredictionet (eval = FALSE)
###################################################
## library(help=predictionet)
###################################################
### code chunk number 9: loadpredictionet (eval = FALSE)
###################################################
## help(expO.colon.ras)
## help(jorissen.colon.ras)
###################################################
### code chunk number 10: helpnetinf (eval = FALSE)
###################################################
## help(netinf)
###################################################
### code chunk number 11: firstnetinf
###################################################
## 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 12: firstnetinf
###################################################
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 13: regrnetopofig
###################################################
## network topology
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 14: edgecoldiff
###################################################
## 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 15: edgecoldiffig
###################################################
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 16: regrnetinfcv
###################################################
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 17: rescv
###################################################
print(str(myres.cv, 1))
###################################################
### code chunk number 18: edgecol
###################################################
## 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 19: edgestabfig0 (eval = FALSE)
###################################################
## 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 20: edgestabfig
###################################################
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: genecolr2
###################################################
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 22: genepredabr2fig0 (eval = FALSE)
###################################################
## 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 23: genepredabr2fig
###################################################
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: genecolmcc
###################################################
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 25: genepredabmccfig0 (eval = FALSE)
###################################################
## 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 26: genepredabmccfig
###################################################
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: validationpred
###################################################
## make the network model predictive
mynet <- net2pred(net=mynet, data=data.ras[ ,goi], method="linear")
## 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 28: genecolr2test
###################################################
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 29: genepredabr2testfig0 (eval = FALSE)
###################################################
## 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 30: genepredabr2testfig
###################################################
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: genecolmcctest
###################################################
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 32: genepredabmcctestfig0 (eval = FALSE)
###################################################
## 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 33: genepredabmcctestfig
###################################################
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: regnetcvexport (eval = FALSE)
###################################################
## netinf2gml(object=myres.cv, file="regrnet_expO_cv")
###################################################
### code chunk number 35: regnetcvpriorsweight
###################################################
## 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 36: regnetcvpriorsweightfig20 (eval = FALSE)
###################################################
## 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 37: regnetcvpriorsweightfig2
###################################################
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: boxplotstabpwcvfig2
###################################################
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 39: boxplotr2pwcvfig2
###################################################
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 40: regnetcvpriorsweightcompr2
###################################################
## 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 41: netinfcv12
###################################################
## 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 42: regnetcvtopo1topo2fig20 (eval = FALSE)
###################################################
## 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 43: regnetcvtopo1topo2fig2
###################################################
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: regnetcvtopo1topo2stabfig20 (eval = FALSE)
###################################################
## 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 45: regnetcvtopo1topo2stabfig2
###################################################
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: regnetcvtopo1topo2predfig2
###################################################
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 47: sessioninfo
###################################################
toLatex(sessionInfo(), locale=FALSE)
Try the predictionet package in your browser
Any scripts or data that you put into this service are public.
predictionet documentation built on Nov. 8, 2020, 7:48 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### R code from vignette source 'predictionet.Rnw'
###################################################
### code chunk number 1: setup
###################################################
## initialize seed
set.seed(123)
###################################################
### code chunk number 2: loadlib
###################################################
library(predictionet)
###################################################
### code chunk number 3: preparepriors
###################################################
## RAS-related genes
genes.ras <- colnames(data.ras)
## read priors generated by the Predictive Networks web application
pn.priors <- read.csv(system.file(file.path("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 4: priorsrandom
###################################################
pn.priors <- pn.priors[sample(1:nrow(pn.priors)), ]
###################################################
### code chunk number 5: tablepriors
###################################################
print(head(pn.priors))
###################################################
### code chunk number 6: preparepriorscount
###################################################
## 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 7: priorscounttab
###################################################
print(table(pn.priors.counts))
###################################################
### code chunk number 8: foopredictionet (eval = FALSE)
###################################################
## library(help=predictionet)
###################################################
### code chunk number 9: loadpredictionet (eval = FALSE)
###################################################
## help(expO.colon.ras)
## help(jorissen.colon.ras)
###################################################
### code chunk number 10: helpnetinf (eval = FALSE)
###################################################
## help(netinf)
###################################################
### code chunk number 11: firstnetinf
###################################################
## 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 12: firstnetinf
###################################################
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 13: regrnetopofig
###################################################
## network topology
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 14: edgecoldiff
###################################################
## 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 15: edgecoldiffig
###################################################
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 16: regrnetinfcv
###################################################
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 17: rescv
###################################################
print(str(myres.cv, 1))
###################################################
### code chunk number 18: edgecol
###################################################
## 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 19: edgestabfig0 (eval = FALSE)
###################################################
## 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 20: edgestabfig
###################################################
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: genecolr2
###################################################
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 22: genepredabr2fig0 (eval = FALSE)
###################################################
## 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 23: genepredabr2fig
###################################################
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: genecolmcc
###################################################
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 25: genepredabmccfig0 (eval = FALSE)
###################################################
## 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 26: genepredabmccfig
###################################################
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: validationpred
###################################################
## make the network model predictive
mynet <- net2pred(net=mynet, data=data.ras[ ,goi], method="linear")
## 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 28: genecolr2test
###################################################
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 29: genepredabr2testfig0 (eval = FALSE)
###################################################
## 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 30: genepredabr2testfig
###################################################
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: genecolmcctest
###################################################
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 32: genepredabmcctestfig0 (eval = FALSE)
###################################################
## 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 33: genepredabmcctestfig
###################################################
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: regnetcvexport (eval = FALSE)
###################################################
## netinf2gml(object=myres.cv, file="regrnet_expO_cv")
###################################################
### code chunk number 35: regnetcvpriorsweight
###################################################
## 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 36: regnetcvpriorsweightfig20 (eval = FALSE)
###################################################
## 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 37: regnetcvpriorsweightfig2
###################################################
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: boxplotstabpwcvfig2
###################################################
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 39: boxplotr2pwcvfig2
###################################################
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 40: regnetcvpriorsweightcompr2
###################################################
## 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 41: netinfcv12
###################################################
## 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 42: regnetcvtopo1topo2fig20 (eval = FALSE)
###################################################
## 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 43: regnetcvtopo1topo2fig2
###################################################
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: regnetcvtopo1topo2stabfig20 (eval = FALSE)
###################################################
## 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 45: regnetcvtopo1topo2stabfig2
###################################################
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: regnetcvtopo1topo2predfig2
###################################################
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 47: sessioninfo
###################################################
toLatex(sessionInfo(), locale=FALSE)
Try the predictionet package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.