################################################### ### chunk number 1: no.nonsense ################################################### #line 42 "vignettes/nem/inst/doc/nem.Rnw" rm(list=ls()) ################################################### ### chunk number 2: ################################################### #line 95 "vignettes/nem/inst/doc/nem.Rnw" library(nem) data("BoutrosRNAi2002") ################################################### ### chunk number 3: ################################################### #line 122 "vignettes/nem/inst/doc/nem.Rnw" res.disc <- nem.discretize(BoutrosRNAiExpression,neg.control=1:4,pos.control=5:8,cutoff=.7) ################################################### ### chunk number 4: ################################################### #line 126 "vignettes/nem/inst/doc/nem.Rnw" disc <- cbind(res.disc$neg,res.disc$pos,res.disc$dat) e.2fold <- BoutrosRNAiExpression[res.disc$sel,] #--- hierarchisch clustern hc <- hclust(as.dist(hamming.distance(disc[,9:16]))) e.2fold <- e.2fold[hc$order, ] disc <- disc [hc$order, ] ################################################### ### chunk number 5: data_cont ################################################### #line 136 "vignettes/nem/inst/doc/nem.Rnw" #--- CONTINUOUS DATA #pdf("data_cont.pdf",width=5,height=13) par(las=2,mgp=c(5.5,1,0),mar=c(6.7,7,4,1),cex.lab=1.7,cex.main=2) image(x = 1:ncol(e.2fold), y = 1:nrow(e.2fold), z = scale(t(e.2fold)), main= "Original data", xlab= "Experiments", xaxt= "n", ylab= "E-genes", yaxt= "n", col = gray(seq(0,.95,length=10)) ) abline(v=c(4,8,10,12,14)+.5) axis(1,1:ncol(e.2fold),colnames(e.2fold)) axis(2,1:nrow(e.2fold),rownames(e.2fold)) #dev.off() ################################################### ### chunk number 6: data_disc ################################################### #line 156 "vignettes/nem/inst/doc/nem.Rnw" #--- DISCRETE DATA #pdf("data_disc.pdf",width=5,height=13) par(las=2,mgp=c(5.5,1,0),mar=c(6.7,7,4,1),cex.lab=1.7,cex.main=2) image(x = 1:ncol(disc), z = t(disc), main= "Discretized data", xlab= "Experiments", xaxt= "n", ylab= "", yaxt= "n", col = gray(seq(.95,0,length=10)) ) abline(v=c(4,8,10,12,14)+.5) axis(1,1:ncol(e.2fold),colnames(e.2fold)) #dev.off() ################################################### ### chunk number 7: ################################################### #line 190 "vignettes/nem/inst/doc/nem.Rnw" res.disc$para ################################################### ### chunk number 8: ################################################### #line 224 "vignettes/nem/inst/doc/nem.Rnw" hyper = set.default.parameters(unique(colnames(res.disc$dat)), para=res.disc$para) result <- nem(res.disc$dat,inference="search",control=hyper, verbose=FALSE) result ################################################### ### chunk number 9: ################################################### #line 233 "vignettes/nem/inst/doc/nem.Rnw" plot.nem(result,what="graph") plot.nem(result,what="mLL") plot.nem(result,what="pos") ################################################### ### chunk number 10: ################################################### #line 239 "vignettes/nem/inst/doc/nem.Rnw" Sgenes <- unique(colnames(res.disc$dat)) models <- enumerate.models(Sgenes) best5 <- -sort(-unique(result$mLL))[1:5] col<-c("yellow","yellow","green","blue") names(col) = Sgenes for (i in 1:5) { graph <- as(models[[which(result$mLL == best5[i])[1]]]-diag(4),"graphNEL") pdf(file=paste("topo",i,".pdf",sep="")) par(cex.main=5) plot(graph, nodeAttrs=list(fillcolor=col), main=paste("-",i, "-")) dev.off() } ################################################### ### chunk number 11: scores1 ################################################### #line 268 "vignettes/nem/inst/doc/nem.Rnw" plot.nem(result,what="mLL") ################################################### ### chunk number 12: pos1 ################################################### #line 275 "vignettes/nem/inst/doc/nem.Rnw" plot.nem(result,what="pos") ################################################### ### chunk number 13: ################################################### #line 293 "vignettes/nem/inst/doc/nem.Rnw" hyper$type="FULLmLL" hyper$hyperpara=c(1,9,9,1) result2 <- nem(res.disc$dat,inference="search",control=hyper,verbose=FALSE) result2 ################################################### ### chunk number 14: ################################################### #line 300 "vignettes/nem/inst/doc/nem.Rnw" best5 <- -sort(-unique(result2$mLL))[1:5] for (i in 1:5) { graph <- as(models[[which(result2$mLL == best5[i])[1]]]-diag(4),"graphNEL") pdf(file=paste("topo2",i,".pdf",sep="")) par(cex.main=5) plot(graph, nodeAttrs=list(fillcolor=col), main=paste("-",i, "-")) dev.off() } ################################################### ### chunk number 15: scores2 ################################################### #line 325 "vignettes/nem/inst/doc/nem.Rnw" plot.nem(result2,what="mLL") ################################################### ### chunk number 16: pos2 ################################################### #line 332 "vignettes/nem/inst/doc/nem.Rnw" plot.nem(result2,what="pos") ################################################### ### chunk number 17: ################################################### #line 351 "vignettes/nem/inst/doc/nem.Rnw" resultPairs <- nem(res.disc$dat,inference="pairwise",control=hyper, verbose=FALSE) resultPairs ################################################### ### chunk number 18: graph3 ################################################### #line 356 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resultPairs$graph) plot.nem(resultPairs, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 19: ################################################### #line 370 "vignettes/nem/inst/doc/nem.Rnw" resultTriples <- nem(res.disc$dat,inference="triples",control=hyper, verbose=FALSE) resultTriples ################################################### ### chunk number 20: graph4 ################################################### #line 375 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resultTriples$graph) plot.nem(resultTriples, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 21: ################################################### #line 387 "vignettes/nem/inst/doc/nem.Rnw" resultGreedy <- nem(res.disc$dat,inference="nem.greedy",control=hyper, verbose=FALSE) resultGreedy ################################################### ### chunk number 22: graph44 ################################################### #line 392 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resultGreedy$graph) plot.nem(resultGreedy, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 23: eval=FALSE ################################################### ## #line 424 "vignettes/nem/inst/doc/nem.Rnw" ## resultMN <- nem(res.disc$dat,inference="ModuleNetwork",control=hyper, verbose=FALSE) # this will do exactly the same as the exhaustive search ################################################### ### chunk number 24: ################################################### #line 433 "vignettes/nem/inst/doc/nem.Rnw" resGreedyMAP <- nem(BoutrosRNAiLods, inference="nem.greedyMAP", control=hyper, verbose=FALSE) resGreedyMAP ################################################### ### chunk number 25: graph55 ################################################### #line 438 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resGreedyMAP$graph) plot.nem(resGreedyMAP, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 26: ################################################### #line 460 "vignettes/nem/inst/doc/nem.Rnw" hyper$Pm = diag(4) hyper$lambda = 10 resultRegularization <- nem(res.disc$dat, inference="search", control=hyper, verbose=FALSE) resultRegularization ################################################### ### chunk number 27: graph6 ################################################### #line 467 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resultRegularization$graph) plot.nem(resultRegularization, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 28: ################################################### #line 490 "vignettes/nem/inst/doc/nem.Rnw" resultModsel <- nemModelSelection(c(0.01,1,100),res.disc$dat, control=hyper,verbose=FALSE) ################################################### ### chunk number 29: graph7 ################################################### #line 496 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resultModsel$graph) plot.nem(resultModsel, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 30: ################################################### #line 507 "vignettes/nem/inst/doc/nem.Rnw" hyper$lambda=0 # set regularization parameter to 0 hyper$Pm # this is our prior graph structure resultBayes <- nem(res.disc$dat, control=hyper, verbose=FALSE) # now we use Bayesian model selection to incorporate the prior information resultBayes ################################################### ### chunk number 31: graph77 ################################################### #line 514 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resultBayes$graph) plot.nem(resultBayes, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 32: eval=FALSE ################################################### ## #line 531 "vignettes/nem/inst/doc/nem.Rnw" ## logdensities = getDensityMatrix(myPValueMatrix,dirname="DiagnosticPlots") ## nem(logdensities[res.disc$sel,],type="CONTmLLBayes",inference="search") ## nem(logdensities[res.disc$sel,],type="CONTmLLMAP",inference="search") ## ## preprocessed <- nem.cont.preprocess(BoutrosRNAiExpression,neg.control=1:4,pos.control=5:8) ## nem(preprocessed$prob.influenced,type="CONTmLL",inference="search") ################################################### ### chunk number 33: eval=FALSE ################################################### ## #line 551 "vignettes/nem/inst/doc/nem.Rnw" ## mydat = filterEGenes(Porig, logdensities) # a-priori filtering of E-genes ## ## hyper$selEGenes = TRUE ## hyper$type = "CONTmLLBayes" ## resAuto = nem(mydat,control=hyper) # use filtered data to estimate a network; perform automated subset selection of E-genes ################################################### ### chunk number 34: eval=FALSE ################################################### ## #line 570 "vignettes/nem/inst/doc/nem.Rnw" ## significance=nem.calcSignificance(disc$dat,res) # assess statistical significance ## bootres=nem.bootstrap(res.disc$dat) # bootstrapping on E-genes ## jackres=nem.jackknife(res.disc$dat) # jackknife on S-genes ## consens=nem.consensus(res.disc$dat) # bootstrap & jackknife on S-genes ################################################### ### chunk number 35: ################################################### #line 581 "vignettes/nem/inst/doc/nem.Rnw" hyper$mode="binary_Bayesian" resultBN = nem(res.disc$dat, inference="BN.greedy", control=hyper) ################################################### ### chunk number 36: graph8 ################################################### #line 586 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resultBN$graph) plot.nem(resultBN, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 37: eval=FALSE ################################################### ## #line 606 "vignettes/nem/inst/doc/nem.Rnw" ## data(SahinRNAi2008) ## control = set.default.parameters(setdiff(colnames(dat.normalized),"time"), map=map.int2node, type="depn",debug=FALSE) # set mapping of interventions to perturbed nodes ## net = nem(dat.normalized, control=control) # greedy hillclimber to find most probable network structure ################################################### ### chunk number 38: ################################################### #line 617 "vignettes/nem/inst/doc/nem.Rnw" resEdgeInf = infer.edge.type(result, BoutrosRNAiLogFC) plot.nem(resEdgeInf, SCC=FALSE) ################################################### ### chunk number 39: graph100 ################################################### #line 623 "vignettes/nem/inst/doc/nem.Rnw" col<-c("yellow","yellow","green","blue") names(col) = nodes(resEdgeInf$graph) plot.nem(resEdgeInf, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 40: ################################################### #line 634 "vignettes/nem/inst/doc/nem.Rnw" plotEffects(res.disc$dat,result) ################################################### ### chunk number 41: plot_effects ################################################### #line 639 "vignettes/nem/inst/doc/nem.Rnw" plotEffects(res.disc$dat,result) ################################################### ### chunk number 42: ################################################### #line 655 "vignettes/nem/inst/doc/nem.Rnw" result.scc <- SCCgraph(result$graph,name=TRUE) plot(result.scc$graph) ################################################### ### chunk number 43: scc ################################################### #line 661 "vignettes/nem/inst/doc/nem.Rnw" # col2<-c("yellow","blue","green") # names(col2) = nodes(result.scc$graph) plot(result.scc$graph)#,nodeAttrs=list(fillcolor=col2)) ################################################### ### chunk number 44: ################################################### #line 687 "vignettes/nem/inst/doc/nem.Rnw" toLatex(sessionInfo())