### R code from vignette source 'GlobalAncova.rnw'

###################################################
### code chunk number 1: GlobalAncova.rnw:39-43
###################################################
#library(Biobase)
library(GlobalAncova)
library(globaltest)
sI <- sessionInfo()


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

data(Golub_Merge)
golubX <- justvsn(Golub_Merge)


###################################################
### code chunk number 3: GlobalAncova.rnw:184-185
###################################################
options(width=70)


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


###################################################
### code 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)


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


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


###################################################
### code chunk number 8: GlobalAncova.rnw:262-263
###################################################
gt.all


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


###################################################
### code chunk number 10: GlobalAncova.rnw:279-280
###################################################
cellcycle <- kegg[["04110"]]


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


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


###################################################
### code chunk number 13: gt.cellcycle
###################################################
gt.cc <- gt(ALL.AML, golubX, subsets=cellcycle)
gt.cc


###################################################
### code chunk number 14: GlobalAncova.rnw:307-308
###################################################
gt.cc


###################################################
### code 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


###################################################
### code chunk number 16: gt.adjust
###################################################
gt.cc.BMPB <- gt(ALL.AML ~ BM.PB, golubX, subsets=cellcycle)
gt.cc.BMPB


###################################################
### code chunk number 17: GlobalAncova.rnw:339-340
###################################################
gt.cc.BMPB


###################################################
### code chunk number 18: GlobalAncova.rnw:350-353
###################################################
data(vantVeer)
data(phenodata)	
data(pathways)	


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


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


###################################################
### code 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


###################################################
### code 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


###################################################
### code 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


###################################################
### code 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


###################################################
### code 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


###################################################
### code chunk number 26: gt.pathways
###################################################
gt.options(transpose = TRUE)
gt.pw <- gt(metastases, vantVeer, subsets=pathways[1:4])
gt.pw


###################################################
### code chunk number 27: GlobalAncova.rnw:503-509
###################################################
ga.pw.raw <- ga.pw[ ,"p.perm"] 
ga.pw.adj <- p.adjust(ga.pw.raw, "holm")
ga.result <- data.frame(rawp=ga.pw.raw, Holm=ga.pw.adj)
ga.result
gt.result <- p.adjust(gt.pw)
gt.result


###################################################
### code 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="")
}


###################################################
### code 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")


###################################################
### code 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


###################################################
### code chunk number 31: gago
###################################################
library(GO.db)
descendants <- get("GO:0007049", GOBPOFFSPRING)
gago <- GAGO(xx=exprs(golubX), formula.full=~ALL.AML, formula.red=~1, model.dat=pData(golubX),
      id=c("GO:0007049", descendants), annotation="hu6800", ontology="BP",
      multtest="focuslevel")


###################################################
### code chunk number 32: GlobalAncova.rnw:617-618
###################################################
head(gago)


###################################################
### code chunk number 33: gakegg
###################################################
GAKEGG(xx=exprs(golubX), formula.full=~ALL.AML, formula.red=~1, model.dat=pData(golubX),
      id=c("04110", "04210"), annotation="hu6800", multtest="BH")


###################################################
### code chunk number 34: gabroad (eval = FALSE)
###################################################
## broad <- getBroadSets("your/path/to/msigdb_v.2.5.xml")
## GABroad(xx=exprs(golubX), formula.full=~ALL.AML, formula.red=~1, model.dat=pData(golubX),
##         id="chr5q33", collection=broad, annotation="hu6800")


###################################################
### code chunk number 35: gabroad2 (eval = FALSE)
###################################################
## GABroad(xx=exprs(golubX), formula.full=~ALL.AML, formula.red=~1, model.dat=pData(golubX),
##         category=c("c1", "c3"), collection=broad, annotation="hu6800", multtest="holm")


###################################################
### code chunk number 36: geneplot
###################################################
Plot.genes(xx=vantVeer, group=metastases, test.genes=p53)
gt.vV <- gt(metastases, vantVeer, subsets=p53)
features(gt.vV, what="w")


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


###################################################
### code chunk number 38: gt_genes
###################################################
features(gt.vV, what="w")


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


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


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


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


###################################################
### code chunk number 43: gt_subjects
###################################################
subjects(gt.vV, what="w")


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


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