# May 1, 2008
# Multiple Testing
# Genome 560

#########################################
# Load libraries for multiple test corrections
library(multtest)
library(qvalue)
#########################################

#########################################
# Load a gene expression data set of 3051 genes and 38 tumor mRNA samples
# 27 acute lymphoblastic leukemia and 11 acute myeloid leukemia
data(golub)
#########################################

#########################################
# use multtest to calculate permutation p-values for the first 1000 genes
smallgd < -golub[1:100,] 
classlabel <- golub.cl
res1 <- mt.maxT(smallgd,classlabel) 
rawp <- res1$rawp[order(res1$index)] 
rawp.500 <- rawp[1:500]
#########################################

#########################################
# now let's use multtest to get adjusted p-values
procs <- c("Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY") 
adj.all <- mt.rawp2adjp(rawp, procs) 
adj.500 <- mt.rawp2adjp(rawp.500, procs) 
#########################################

#########################################
# compare number of bonferroni corrected p-values significant at p < 0.05 between the first 500 and full 1000 genes
length(which(adj.500$adjp[,2] <= 0.05))
length(which(adj.all$adjp[,2] <= 0.05))
#########################################

#########################################
# compare results by plotting raw versus adjusted p-values
plot(adj.all$adjp[,1], adj.all$adjp[,2], pch = 19)  #Bonferroni
points(adj.all$adjp[,1], adj.all$adjp[,3], col = "red", pch = 19) #Holm
points(adj.all$adjp[,1], adj.all$adjp[,7], col = "green", pch = 19) #BH
#########################################

#########################################
# qvalue
#########################################
qobj <- qvalue(rawp)
summary(qobj, cuts=c(0.01, 0.05))
qplot(qobj)
qobj$qvalues
qwrite(qobj, filename="myresults.txt")
#########################################

#########################################
# Compare pFDR (qvalues) with BH FDRs
qval <- qobj$qvalues
pval <- qobj$pvalues
length(which(adj.all$adjp[,7] <= 0.05))
length(which(qval <= 0.05))

# get bh fdrs for original order of rawp and plot it versus qvalue
bh <- adj.all$adjp[order(adj.all$index),7]
plot(bh, qval, pch = 19)
abline(h=0.5)
abline(v=0.05)
#########################################