### R code from vignette source 'vignettes/predictionet/inst/doc/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)