###################################################
### chunk number 1: 
###################################################
#library(Biobase)
library(GlobalAncova)
library(globaltest)
sI <- sessionInfo()


###################################################
### chunk number 2: initialize
###################################################
library(GlobalAncova)
library(globaltest)
library(golubEsets)
library(hu6800.db)
library(vsn)
library(multtest)

data(Golub_Merge)
golubX <- justvsn(Golub_Merge)


###################################################
### chunk number 3: 
###################################################
options(width=70)


###################################################
### chunk number 4: all
###################################################
gr <- as.numeric(golubX$ALL.AML=="ALL")
ga.all <- GlobalAncova(xx=exprs(golubX), group=gr, covars=NULL, perm=100)


###################################################
### chunk number 5: all2 eval=FALSE
###################################################
## GlobalAncova(xx=exprs(golubX), formula.full=~ALL.AML, formula.red=~1, model.dat=pData(golubX), perm=100)
## GlobalAncova(xx=exprs(golubX), formula.full=~ALL.AML, test.terms="ALL.AMLAML", model.dat=pData(golubX), perm=100)


###################################################
### chunk number 6: all.display
###################################################
ga.all 


###################################################
### chunk number 7: gt.all
###################################################
gt.all <- globaltest(golubX, "ALL.AML")
gt.all


###################################################
### chunk number 8: 
###################################################
gt.all


###################################################
### chunk number 9: loadKEGG
###################################################
kegg <- as.list(hu6800PATH2PROBE)


###################################################
### chunk number 10: 
###################################################
cellcycle <- kegg[["04110"]]


###################################################
### chunk number 11: cc eval=FALSE
###################################################
## cellcycle <- kegg[["04110"]]


###################################################
### chunk number 12: cellcycle
###################################################
ga.cc <- GlobalAncova(xx=exprs(golubX), group=gr, test.genes=cellcycle, method="both", perm=1000) 
ga.cc


###################################################
### chunk number 13: gt.cellcycle
###################################################
gt.cc <- globaltest(X=golubX, Y="ALL.AML", genesets=cellcycle)
gt.cc


###################################################
### chunk number 14: 
###################################################
gt.cc


###################################################
### chunk number 15: adjust
###################################################
ga.cc.BMPB <- GlobalAncova(xx=exprs(golubX), group=gr, covars=golubX$BM.PB, 
test.genes=cellcycle, method="both", perm=1000) 
ga.cc.BMPB


###################################################
### chunk number 16: gt.adjust
###################################################
gt.cc.BMPB <- globaltest(X=golubX, Y=ALL.AML ~ BM.PB, genesets=cellcycle)
gt.cc.BMPB


###################################################
### chunk number 17: 
###################################################
gt.cc.BMPB


###################################################
### chunk number 18: 
###################################################
data(vantVeer)
data(phenodata)	
data(pathways)	


###################################################
### chunk number 19: vV
###################################################
data(vantVeer)
data(phenodata)	
data(pathways)
metastases <- phenodata$metastases
p53 <- pathways$p53_signalling	


###################################################
### chunk number 20: p53test
###################################################
vV.1 <- GlobalAncova(xx=vantVeer, group=metastases, test.genes=p53, method="both", perm=1000)
vV.1


###################################################
### chunk number 21: grade
###################################################
vV.3 <- GlobalAncova(xx=vantVeer, formula.full=~grade, formula.red=~1, model.dat=phenodata, test.genes=p53, method="both", perm=1000)
vV.3


###################################################
### chunk number 22: interact
###################################################
signature.gene <- "AL137718"
model <- data.frame(phenodata, signature.gene=vantVeer[signature.gene, ])
vV.4 <- GlobalAncova(xx=vantVeer, formula.full=~signature.gene + ERstatus, formula.red=~ERstatus, model.dat=model, test.genes=p53, method="both", perm=1000)
vV.4


###################################################
### chunk number 23: coexpr
###################################################
vV.5 <- GlobalAncova(xx=vantVeer, formula.full=~metastases + signature.gene + ERstatus, formula.red=~metastases + ERstatus, model.dat=model, test.genes=p53, method="both", perm=1000)
vV.5


###################################################
### chunk number 24: diffcoexpr
###################################################
vV.6 <- GlobalAncova(xx=vantVeer, formula.full=~metastases * signature.gene + ERstatus, formula.red=~metastases + signature.gene + ERstatus, model.dat=model, test.genes=p53, method="both", perm=1000)
vV.6


###################################################
### chunk number 25: pathways
###################################################
metastases <- phenodata$metastases
ga.pw <- GlobalAncova(xx=vantVeer, group=metastases, test.genes=pathways[1:4], method="both", perm=1000) 
ga.pw


###################################################
### chunk number 26: gt.pathways
###################################################
gt.pw <- globaltest(X=vantVeer, Y=metastases, genesets=pathways[1:4])
gt.pw


###################################################
### chunk number 27: 
###################################################
ga.pw.raw <- ga.pw[ ,"p.perm"]
ga.pw.adj <- mt.rawp2adjp(ga.pw.raw)
ga.result <- ga.pw.adj$adjp[order(ga.pw.adj$index), c("rawp", "Holm")]
rownames(ga.result) <- names(pathways)[1:4]
ga.result             
gt.pw.raw <- p.value(gt.pw)   
gt.pw.adj <- mt.rawp2adjp(gt.pw.raw)
gt.result <- gt.pw.adj$adjp[order(gt.pw.adj$index), c("rawp", "Holm")]
rownames(gt.result) <- names(pathways)[1:4]             
gt.result


###################################################
### chunk number 28: hypoth
###################################################
if(require(Rgraphviz))
{
res <- GlobalAncova:::Hnull.family(pathways[1:4])
graph <- new("graphNEL", nodes=names(res), edgemode="directed")
graph <- addEdge(from=c(rep(names(res)[15],4),rep(names(res)[10],3),rep(names(res)[11],3),
  rep(names(res)[13],3),rep(names(res)[14],3),rep(names(res)[5],2),rep(names(res)[6],2),
  rep(names(res)[7],2),rep(names(res)[8],2),rep(names(res)[9],2),rep(names(res)[12],2)), 
  to=c(names(res)[10],names(res)[11],names(res)[13],names(res)[14],names(res)[5],names(res)[6],
  names(res)[8],names(res)[5],names(res)[7],names(res)[9],names(res)[6],names(res)[7],names(res)[12],
  names(res)[8],names(res)[9],names(res)[12],names(res)[1],names(res)[2],names(res)[1],names(res)[3],
  names(res)[1],names(res)[4],names(res)[2],names(res)[3],names(res)[2],names(res)[4],names(res)[3],
  names(res)[4]), graph, weights=rep(1, 28))

att <- list()
lab <- c("1","2","3","4","1-2","1-3","1-4","2-3","2-4","1-2-3","1-2-4","3-4","1-3-4","2-3-4","1-2-3-4")
names(lab) <- names(res)
att$label <- lab

plot(graph, nodeAttrs=att)
} else {
plot(1, 1, type="n", main="This plot requires Rgraphviz", 
     axes=FALSE, xlab="", ylab="")
}


###################################################
### chunk number 29: closed.test
###################################################
ga.closed <- GlobalAncova.closed(xx=vantVeer, group=metastases, 
test.genes=pathways[1:4], previous.test=ga.pw, level=0.05, method="approx")


###################################################
### chunk number 30: names.closed.test
###################################################
names(ga.closed)
rownames(ga.closed$test.results[[1]])
rownames(ga.closed$test.results[[1]]) <- NULL
ga.closed$test.results[1]
ga.closed$significant
ga.closed$not.significant


###################################################
### chunk number 31: makeGOstructure
###################################################
bp <- makeGOstructure(golubX, "hu6800")


###################################################
### chunk number 32: getFocus
###################################################
focusBP <- getFocus(bp, maxatoms=5)
str(focusBP)


###################################################
### chunk number 33: gtGO
###################################################
go30 <- GAGO(exprs(golubX), group=pData(golubX)$ALL.AML, focus = focusBP, GO = bp, stopafter = 30)


###################################################
### chunk number 34: 
###################################################
str(go30)


###################################################
### chunk number 35: geneplot
###################################################
Plot.genes(xx=vantVeer, group=metastases, test.genes=p53)
gt.vV <- globaltest(X=vantVeer, Y=metastases, genesets=p53)
geneplot(gt.vV)


###################################################
### chunk number 36: genes
###################################################
Plot.genes(xx=vantVeer, group=metastases, test.genes=p53)


###################################################
### chunk number 37: gt_genes
###################################################
geneplot(gt.vV)


###################################################
### chunk number 38: geneplot
###################################################
Plot.genes(xx=vantVeer, formula.full=~metastases, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="ERstatus")


###################################################
### chunk number 39: genes2
###################################################
Plot.genes(xx=vantVeer, formula.full=~metastases, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="ERstatus")


###################################################
### chunk number 40: subjectsplot
###################################################
#colnames(exprs(golubX)) <- pData(golubX)[ ,1]
Plot.subjects(xx=vantVeer, group=metastases, test.genes=p53, legendpos="bottomright")
sampleplot(gt.vV)


###################################################
### chunk number 41: subjects
###################################################
Plot.subjects(xx=vantVeer, group=metastases, test.genes=p53, legendpos="bottomright")


###################################################
### chunk number 42: gt_subjects
###################################################
sampleplot(gt.vV)


###################################################
### chunk number 43: subjectsplot2
###################################################
Plot.subjects(xx=vantVeer, formula.full=~grade, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="grade", legendpos="topleft")


###################################################
### chunk number 44: subjects2
###################################################
Plot.subjects(xx=vantVeer, formula.full=~grade, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="grade", legendpos="topleft")