################################################### ### chunk number 1: no.nonsense ################################################### rm(list=ls()) ################################################### ### chunk number 2: ################################################### library(nem) data("BoutrosRNAi2002") ################################################### ### chunk number 3: ################################################### res.disc <- nem.discretize(BoutrosRNAiExpression,neg.control=1:4,pos.control=5:8,cutoff=.7) ################################################### ### chunk number 4: ################################################### 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 ################################################### #--- 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 ################################################### #--- 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: ################################################### res.disc$para ################################################### ### chunk number 8: ################################################### 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: ################################################### plot.nem(result,what="graph") plot.nem(result,what="mLL") plot.nem(result,what="pos") ################################################### ### chunk number 10: ################################################### 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 ################################################### plot.nem(result,what="mLL") ################################################### ### chunk number 12: pos1 ################################################### plot.nem(result,what="pos") ################################################### ### chunk number 13: ################################################### 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: ################################################### 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 ################################################### plot.nem(result2,what="mLL") ################################################### ### chunk number 16: pos2 ################################################### plot.nem(result2,what="pos") ################################################### ### chunk number 17: ################################################### resultPairs <- nem(res.disc$dat,inference="pairwise",control=hyper, verbose=FALSE) resultPairs ################################################### ### chunk number 18: graph3 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultPairs$graph) plot.nem(resultPairs, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 19: ################################################### resultTriples <- nem(res.disc$dat,inference="triples",control=hyper, verbose=FALSE) resultTriples ################################################### ### chunk number 20: graph4 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultTriples$graph) plot.nem(resultTriples, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 21: ################################################### resultGreedy <- nem(res.disc$dat,inference="nem.greedy",control=hyper, verbose=FALSE) resultGreedy ################################################### ### chunk number 22: graph44 ################################################### 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 ################################################### ## resultMN <- nem(res.disc$dat,inference="ModuleNetwork",control=hyper, verbose=FALSE) # this will do exactly the same as the exhaustive search ################################################### ### chunk number 24: ################################################### resGreedyMAP <- nem(BoutrosRNAiLods, inference="nem.greedyMAP", control=hyper, verbose=FALSE) resGreedyMAP ################################################### ### chunk number 25: graph55 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resGreedyMAP$graph) plot.nem(resGreedyMAP, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 26: ################################################### hyper$Pm = diag(4) hyper$lambda = 10 resultRegularization <- nem(res.disc$dat, inference="search", control=hyper, verbose=FALSE) resultRegularization ################################################### ### chunk number 27: graph6 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultRegularization$graph) plot.nem(resultRegularization, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 28: ################################################### resultModsel <- nemModelSelection(c(0.01,1,100),res.disc$dat, control=hyper,verbose=FALSE) ################################################### ### chunk number 29: graph7 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultModsel$graph) plot.nem(resultModsel, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 30: ################################################### 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 ################################################### 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 ################################################### ## 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 ################################################### ## 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 ################################################### ## 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: ################################################### hyper$mode="binary_Bayesian" resultBN = nem(res.disc$dat, inference="BN.greedy", control=hyper) ################################################### ### chunk number 36: graph8 ################################################### 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 ################################################### ## control = set.default.parameters(setdiff(colnames(mydata),"time"), map=mymap, trans.close=FALSE, type="gnem",debug=FALSE, Pm=diag(10)) # set mapping of interventions to perturbed nodes; do not require transitive closed graphs; employ backward elimination for greedy hillclimbing or module networks; require sparsity of the network ## ## net = nem(dat, control=control) # greedy hillclimber with final backward elimination step ################################################### ### chunk number 38: ################################################### resEdgeInf = infer.edge.type(result, BoutrosRNAiLogFC) plot.nem(resEdgeInf, SCC=FALSE) ################################################### ### chunk number 39: graph100 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resEdgeInf$graph) plot.nem(resEdgeInf, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### chunk number 40: ################################################### plotEffects(res.disc$dat,result) ################################################### ### chunk number 41: plot_effects ################################################### plotEffects(res.disc$dat,result) ################################################### ### chunk number 42: ################################################### result.scc <- SCCgraph(result$graph,name=TRUE) plot(result.scc$graph) ################################################### ### chunk number 43: scc ################################################### # col2<-c("yellow","blue","green") # names(col2) = nodes(result.scc$graph) plot(result.scc$graph)#,nodeAttrs=list(fillcolor=col2)) ################################################### ### chunk number 44: ################################################### toLatex(sessionInfo())