###################################################
### chunk number 1: loadPacks
###################################################
library("RbcBook1")
library("Biobase")
library("multtest")
library("genefilter")


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


###################################################
### chunk number 3: ALL
###################################################
library("ALL")
library("hgu95av2")
data(ALL)


###################################################
### chunk number 4: genefilter
###################################################
ffun <- filterfun(pOverA(p=0.2, A=100), cv(a=0.7, b=10))
filt <- genefilter(2^exprs(ALL), ffun)
filtALL <- ALL[filt,]


###################################################
### chunk number 5: ALLDetails
###################################################
filtX <- exprs(filtALL)
pheno <- pData(filtALL)


###################################################
### chunk number 6: BT
###################################################
table(pData(ALL)$BT)
Bcell<-rep(0,length(pData(ALL)$BT))
Bcell[grep("B",as.character(pData(ALL)$BT))]<-1


###################################################
### chunk number 7: BTBoot
###################################################
seed <- 99
BT.boot <- cache("BT.boot", MTP(X=filtX, Y=Bcell, 
alternative="greater", B=100, method="sd.minP", seed=seed))


###################################################
### chunk number 8: BTOut
###################################################
summary(BT.boot)


###################################################
### chunk number 9: BTGenes
###################################################
BT.diff <- BT.boot@adjp<=0.05
BT.AffyID <- geneNames(filtALL)[BT.diff]
mget(BT.AffyID[1:2], env=hgu95av2GENENAME)


###################################################
### chunk number 10: BTPlot
###################################################
par(mfrow=c(2,2))
plot(BT.boot)


###################################################
### chunk number 11: BTtppfp
###################################################
q <- c(0.05,0.1,0.25)
BT.tppfp <- fwer2tppfp(adjp=BT.boot@adjp, q=q)
comp.tppfp <- cbind(BT.boot@adjp, BT.tppfp)
mtps <- c("FWER",paste("TPPFP(",q,")", sep=""))
mt.plot(adjp=comp.tppfp, teststat=BT.boot@statistic, proc=mtps, 
leg=c(0.1,430), col=1:4, lty=1:4, lwd=3)
title("Comparison of TPPFP(q)-controlling AMTPs\n based on SD minP MTP")


###################################################
### chunk number 12: BTfdr
###################################################
BT.fdr <- fwer2fdr(adjp=BT.boot@adjp, method="both")$adjp
BT.marg.fdr <- mt.rawp2adjp(rawp=BT.boot@rawp, proc=c("BY","BH"))
comp.fdr <- cbind(BT.fdr, BT.marg.fdr$adjp[order(BT.marg.fdr$index),-1])
mtps <- c("AMTP Cons", "AMTP Rest", "BY", "BH")
mt.plot(adjp=comp.fdr, teststat=BT.boot@statistic, proc=mtps, 
leg=c(0.1,430), col=c(2,2,3,3), lty=rep(1:2,2), lwd=3)
title("Comparison of FDR-controlling MTPs")


###################################################
### chunk number 13: mbBoot
###################################################
BX<-filtX[,Bcell==1]
Bpheno<-pheno[Bcell==1,]
mb <- as.character(Bpheno$mol.biol)
table(mb)
other <- c("E2A/PBX1", "p15/p16")
mb.boot <- cache("mb.boot", MTP(X=BX[,!(mb%in%other)], Y=mb[!(mb%in%other)], 
test="f", alpha=c(0.01,0.1), B=100, get.cutoff=TRUE, seed=seed))


###################################################
### chunk number 14: mbOut
###################################################
mb.rej<-summary(mb.boot)$rejections


###################################################
### chunk number 15: mbOut
###################################################
mb.rej


###################################################
### chunk number 16: coxphPrep
###################################################
# Patients with original complete remission and who were followed up
cr.ind <- (Bpheno$remission=="CR")
cr.pheno <- Bpheno[cr.ind,]
times <- strptime(cr.pheno$"date last seen", "%m/%d/%Y")-strptime(cr.pheno$date.cr, "%m/%d/%Y")
time.ind <- !is.na(times)
times <- times[time.ind]
# Patients who haven't relapsed are treated as censored
cens <- ((1:length(times))%in%grep("CR", cr.pheno[time.ind,"f.u"]))
# Time to relapse
rel.times <- Surv(times, !cens)
patients <- (1: ncol(BX))[cr.ind][time.ind]
# Prepare data for MTP
relX <- BX[, patients]
relZ <- Bpheno[patients,]


###################################################
### chunk number 17: coxphBoot
###################################################
cox.boot <- cache("cox.boot", MTP(X=relX, Y=rel.times, Z=relZ,
       Z.incl="sex", Z.test=NULL, test="coxph.YvsXZ", B=100, get.cr=TRUE,
       seed=seed))


###################################################
### chunk number 18: coxphOut
###################################################
cox.diff <- cox.boot@adjp<=0.05
sum(cox.diff)
cox.AffyID <- geneNames(filtALL)[cox.diff]
mget(cox.AffyID, env=hgu95av2GENENAME)


###################################################
### chunk number 19: coxphPlot
###################################################
plot(cox.boot, which=5, top=5, sub.caption=NULL)
abline(h=0,col="red")