###################################################
### 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
###################################################
## 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: 
###################################################
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())