### R code from vignette source 'vignettes/nem/inst/doc/nem.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: no.nonsense ################################################### rm(list=ls()) ################################################### ### code chunk number 2: nem.Rnw:94-96 ################################################### library(nem) data("BoutrosRNAi2002") ################################################### ### code chunk number 3: nem.Rnw:121-122 ################################################### res.disc <- nem.discretize(BoutrosRNAiExpression,neg.control=1:4,pos.control=5:8,cutoff=.7) ################################################### ### code chunk number 4: nem.Rnw:125-132 ################################################### 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, ] ################################################### ### code 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() ################################################### ### code 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() ################################################### ### code chunk number 7: nem.Rnw:189-190 ################################################### res.disc$para ################################################### ### code chunk number 8: nem.Rnw:223-226 ################################################### 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 ################################################### ### code chunk number 9: nem.Rnw:232-235 ################################################### plot.nem(result,what="graph") plot.nem(result,what="mLL") plot.nem(result,what="pos") ################################################### ### code chunk number 10: nem.Rnw:238-252 ################################################### 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() } ################################################### ### code chunk number 11: scores1 ################################################### plot.nem(result,what="mLL") ################################################### ### code chunk number 12: pos1 ################################################### plot.nem(result,what="pos") ################################################### ### code chunk number 13: nem.Rnw:292-296 ################################################### hyper$type="FULLmLL" hyper$hyperpara=c(1,9,9,1) result2 <- nem(res.disc$dat,inference="search",control=hyper,verbose=FALSE) result2 ################################################### ### code chunk number 14: nem.Rnw:299-309 ################################################### 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() } ################################################### ### code chunk number 15: scores2 ################################################### plot.nem(result2,what="mLL") ################################################### ### code chunk number 16: pos2 ################################################### plot.nem(result2,what="pos") ################################################### ### code chunk number 17: nem.Rnw:350-352 ################################################### resultPairs <- nem(res.disc$dat,inference="pairwise",control=hyper, verbose=FALSE) resultPairs ################################################### ### code chunk number 18: graph3 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultPairs$graph) plot.nem(resultPairs, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 19: nem.Rnw:369-371 ################################################### resultTriples <- nem(res.disc$dat,inference="triples",control=hyper, verbose=FALSE) resultTriples ################################################### ### code chunk number 20: graph4 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultTriples$graph) plot.nem(resultTriples, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 21: nem.Rnw:386-388 ################################################### resultGreedy <- nem(res.disc$dat,inference="nem.greedy",control=hyper, verbose=FALSE) resultGreedy ################################################### ### code chunk number 22: graph44 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultGreedy$graph) plot.nem(resultGreedy, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 23: nem.Rnw:423-424 (eval = FALSE) ################################################### ## resultMN <- nem(res.disc$dat,inference="ModuleNetwork",control=hyper, verbose=FALSE) # this will do exactly the same as the exhaustive search ################################################### ### code chunk number 24: nem.Rnw:432-434 ################################################### resGreedyMAP <- nem(BoutrosRNAiLods, inference="nem.greedyMAP", control=hyper, verbose=FALSE) resGreedyMAP ################################################### ### code chunk number 25: graph55 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resGreedyMAP$graph) plot.nem(resGreedyMAP, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 26: nem.Rnw:459-463 ################################################### hyper$Pm = diag(4) hyper$lambda = 10 resultRegularization <- nem(res.disc$dat, inference="search", control=hyper, verbose=FALSE) resultRegularization ################################################### ### code chunk number 27: graph6 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultRegularization$graph) plot.nem(resultRegularization, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 28: nem.Rnw:489-490 ################################################### resultModsel <- nemModelSelection(c(0.01,1,100),res.disc$dat, control=hyper,verbose=FALSE) ################################################### ### code chunk number 29: graph7 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultModsel$graph) plot.nem(resultModsel, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 30: nem.Rnw:506-510 ################################################### 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 ################################################### ### code chunk number 31: graph77 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultBayes$graph) plot.nem(resultBayes, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 32: nem.Rnw:530-536 (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") ################################################### ### code chunk number 33: nem.Rnw:550-557 (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 with tuning of the parameter delta ## ## most.relevant = getRelevantEGenes(as(resAuto$graph, "matrix"), mydat)$selected # returns all E-genes with partial log-likelihood > 0. ################################################### ### code chunk number 34: nem.Rnw:571-575 (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 ################################################### ### code chunk number 35: nem.Rnw:582-584 ################################################### hyper$mode="binary_Bayesian" resultBN = nem(res.disc$dat, inference="BN.greedy", control=hyper) ################################################### ### code chunk number 36: graph8 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultBN$graph) plot.nem(resultBN, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 37: nem.Rnw:598-603 ################################################### hyper$mcmc.nburnin = 10 # much to few in practice hyper$mcmc.nsamples = 100 hyper$type = "CONTmLLDens" hyper$Pm = NULL resulteminem = nem(BoutrosRNAiLods, inference="mc.eminem", control=hyper) ################################################### ### code chunk number 38: graph8 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resultBN$graph) plot.nem(resulteminem, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 39: nem.Rnw:618-623 (eval = FALSE) ################################################### ## data("Ivanova2006RNAiTimeSeries") ## dim(dat) ## control = set.default.parameters(dimnames(dat)[[3]], para=c(0.1,0.1)) ## net = nem(dat, inference="dynoNEM", control=control) ## plot.nem(net, SCC=FALSE, plot.probs=TRUE) ################################################### ### code chunk number 40: nem.Rnw:641-644 (eval = FALSE) ################################################### ## 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 ################################################### ### code chunk number 41: nem.Rnw:652-654 ################################################### resEdgeInf = infer.edge.type(result, BoutrosRNAiLogFC) plot.nem(resEdgeInf, SCC=FALSE) ################################################### ### code chunk number 42: graph100 ################################################### col<-c("yellow","yellow","green","blue") names(col) = nodes(resEdgeInf$graph) plot.nem(resEdgeInf, nodeAttrs=list(fillcolor=col), SCC=FALSE) ################################################### ### code chunk number 43: nem.Rnw:668-669 ################################################### plotEffects(res.disc$dat,result) ################################################### ### code chunk number 44: plot_effects ################################################### plotEffects(res.disc$dat,result) ################################################### ### code chunk number 45: nem.Rnw:689-691 ################################################### result.scc <- SCCgraph(result$graph,name=TRUE) plot(result.scc$graph) ################################################### ### code chunk number 46: scc ################################################### # col2<-c("yellow","blue","green") # names(col2) = nodes(result.scc$graph) plot(result.scc$graph)#,nodeAttrs=list(fillcolor=col2)) ################################################### ### code chunk number 47: nem.Rnw:721-722 ################################################### toLatex(sessionInfo())