###################################################
### chunk number 1: loadPacks
###################################################
library(Biobase)
library(annotate)
library(tkWidgets)
library(genefilter)
library(multtest)


###################################################
### chunk number 2: loadALL
###################################################
library(ALL)
library(hgu95av2)
data(ALL)
class(ALL)
slotNames(ALL)
show(ALL)
varLabels(ALL)


###################################################
### chunk number 3: filterALL
###################################################
X<-exprs(ALL)
## First, build the filters
f1 <- pOverA(p=0.2, A=100)
f2 <- cv(a=0.7, b=10)
## Next, assemble them in a filtering function
ff <- filterfun(f1,f2)
## Then, apply the filtering function
geneSub<-genefilter(2^X,ff)
## Obtain a microarray object for the relevant gene subset
subALL<-ALL[geneSub,]
show(subALL)


###################################################
### chunk number 4: filtALL
###################################################
data(filtALL)
show(filtALL)


###################################################
### chunk number 5: funcMTP
###################################################
args(MTP)


###################################################
### chunk number 6: classMTP
###################################################
slotNames("MTP")


###################################################
### chunk number 7: cytoBoot
###################################################
set.seed(999)
cyto.boot<-MTP(X=exprs(filtALL), Y=pData(filtALL)$cyto.normal, alternative="less", B=500, method="sd.minP")


###################################################
### chunk number 8: cytoOut
###################################################
class(cyto.boot)
slotNames(cyto.boot)
print(cyto.boot)
summ<-summary(cyto.boot) 
class(summ)
names(summ)
summ


###################################################
### chunk number 9: cytoGenes
###################################################
sum(cyto.boot@adjp<=0.05)
genes<-geneNames(filtALL)[cyto.boot@adjp<=0.05]
genes


###################################################
### chunk number 10: cytoAnnotate
###################################################
hgu95av2()
mget(genes, envir=hgu95av2GENENAME)
mget(genes, envir=hgu95av2MAP)
mget(genes, envir=hgu95av2PMID)


###################################################
### chunk number 11: cytoPubMed
###################################################
allAbsts <- pm.getabst(genes, "hgu95av2")
pm.titles(allAbsts)
oneAbst<-allAbsts[[1]][[1]]
class(oneAbst)
slotNames(oneAbst)
oneAbst
abstText(oneAbst)
allAbsts2<-NULL
for(l in allAbsts)
  allAbsts2<-c(allAbsts2,l)
pmAbst2HTML(allAbsts2, filename="cytoBoot", title="Chiaretti et al. (2004) ALL dataset: FWER-controlling bootstrap step-down minP MTP, cytogenetic test status", frames=TRUE)


###################################################
### chunk number 12: cytoPlot
###################################################
par(mfrow=c(2,2))
plot(cyto.boot)


###################################################
### chunk number 13: cytoMarg
###################################################
cyto.marg<-mt.rawp2adjp(rawp=cyto.boot@rawp,proc=c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD"))
compFWER<-cbind(cyto.boot@adjp,cyto.marg$adjp[order(cyto.marg$index),-1])


###################################################
### chunk number 14: mtplot
###################################################
par(mfrow=c(1,1))
mt.plot(adjp=compFWER,teststat=cyto.boot@statistic,proc=c("SD minP","Bonferroni","Holm","Hochberg","SidakSS","SidakSD"), leg=c(0.1,400), col=1:6, lty=1:6, lwd=2)
title("Comparison of FWER-controlling marginal MTPs and step-down minP MTP")


###################################################
### chunk number 15: cytoPermOld
###################################################
NAs<-is.na(pData(filtALL)$cyto.normal)
set.seed(999)
cyto.perm.old<-mt.minP(X=exprs(filtALL)[,!NAs], classlabel=pData(filtALL)$cyto.normal[!NAs], side="lower", B=500)
names(cyto.perm.old)
sum(cyto.perm.old$adjp<0.05)


###################################################
### chunk number 16: cytoPermNew
###################################################
set.seed(999)
cyto.perm.new<-MTP(X=exprs(filtALL), Y=pData(filtALL)$cyto.normal, alternative="less", nulldist="perm", B=500, method="sd.minP")
summary(cyto.perm.new)


###################################################
### chunk number 17: cytoWilcoxon
###################################################
set.seed(999)
cyto.wilcox<-MTP(X=exprs(filtALL), Y=pData(filtALL)$cyto.normal, robust=TRUE, alternative="less", B=500, method="sd.minP")
sum(cyto.wilcox@adjp<0.05)
sum(cyto.wilcox@adjp<0.05 & cyto.boot@adjp<0.05)


###################################################
### chunk number 18: cytogfwer
###################################################
k<-c(10,50,100)
cyto.gfwer<-fwer2gfwer(adjp=cyto.boot@adjp, k=k)
comp.gfwer<-cbind(cyto.boot@adjp,cyto.gfwer)
mtps<-paste("gFWER(",c(0,k),")",sep="")
mt.plot(adjp=comp.gfwer, teststat=cyto.boot@statistic, proc=mtps, leg=c(0.1,400),col=1:4,lty=1:4, lwd=2)
title("Comparison of gFWER(k)-controlling AMTPs based on SD minP")


###################################################
### chunk number 19: cytotppfp
###################################################
q<-c(0.05,0.1,0.5)
cyto.tppfp<-fwer2tppfp(adjp=cyto.boot@adjp, q=q)
comp.tppfp<-cbind(cyto.boot@adjp,cyto.tppfp)
mtps<-c("FWER",paste("TPPFP(",q,")",sep=""))
mt.plot(adjp=comp.tppfp, teststat=cyto.boot@statistic, proc=mtps, leg=c(0.1,400), col=1:4, lty=1:4, lwd=2)
title("Comparison of TPPFP(q)-controlling AMTPs based on SD minP")


###################################################
### chunk number 20: cytofdr
###################################################
set.seed(999)
cyto.fdr<-MTP(X=exprs(filtALL), Y=pData(filtALL)$cyto.normal, alternative="less", typeone="fdr", B=500, method="sd.minP")
cyto.fdr.marg<-mt.rawp2adjp(rawp=cyto.fdr@rawp,proc=c("BY","BH"))
comp.fdr<-cbind(cyto.fdr@adjp,cyto.fdr.marg$adjp[order(cyto.fdr.marg$index),-1])
mt.plot(adjp=comp.fdr,teststat=cyto.fdr@statistic,proc=c("AMTP","BY","BH"),leg=c(0.1,400),col=1:3,lty=1:3, lwd=2)
title("Comparison of FDR-controlling MTPs")


###################################################
### chunk number 21: mbBoot
###################################################
set.seed(999)
mb<-as.character(pData(filtALL)$mol.biol)
other<-c("E2A/PBX1","NUP-98","p15/p16")
mb[mb%in%other]<-"other"
mb.boot<-MTP(X=exprs(filtALL), Y=mb, test="f", alpha=c(0.01,0.1), B=100, get.cutoff=TRUE)


###################################################
### chunk number 22: mbOut
###################################################
sum(mb.boot@adjp<=0.01)
sum(mb.boot@statistic>mb.boot@cutoff[,"alpha=0.01"])


###################################################
### chunk number 23: coxphPrep
###################################################
# Which genes are associated with tumor molecular subtype
genes<-mb.boot@adjp<=0.01
# Which patients have complete remission
cr.ind<-pData(filtALL)$remission=="CR"
crPheno<-pData(filtALL)[cr.ind,]
times<-strptime(crPheno$"date last seen","%m/%d/%Y")-strptime(crPheno$date.cr,"%m/%d/%Y")
times.ind<-!is.na(times)
times<-times[times.ind]
# Which patients still haven't relapsed
cens<-((1:length(times))%in%grep("CR",crPheno[times.ind,"f.u"]))
rel.times<-Surv(times,!cens)
patients<-(1: ncol(exprs(filtALL)))[cr.ind][times.ind]
# Prepare data for MTP
X<-exprs(filtALL[genes, patients])
Z<-pData(filtALL[genes, patients])


###################################################
### chunk number 24: coxphBoot
###################################################
set.seed(999)
cox.boot<-MTP(X=X, Y=rel.times, Z=Z, Z.incl="sex", Z.test=1, test="coxph.YvsXZ", B=200, get.cr=TRUE)
summary(cox.boot)
geneNames(filtALL)[cox.boot@adjp<=0.05]


###################################################
### chunk number 25: coxphPlot
###################################################
plot(cox.boot, which = 5)