empirical_FDR.r

## LOADING LIBRARIES FOR PARALLEL PROCESSING
library(doMC)
ncore = multicore:::detectCores()
registerDoMC(cores = ncore)

## COMPUTING THE RANDOM P-VALUE DISTRIBUTION
# How many random sampling
R=100
# Shuffling the sample labels and recomputing the p-value each time
p.rand <- foreach(i = 1:dim(e)[1], .combine=rbind) %dopar% {
  pval <- vector(mode='numeric', length=R)
  for(j in 1:R){
		randomOrder <- sample(1:length(classLabel))
		pval[j] <- t.test(as.numeric(e[i,randomOrder[1:9]]), as.numeric(e[i,randomOrder[10:18]]))$p.value
	}
	pval
}

## COMPUTING Q-VALUES
qval <- foreach(i = 1:length(p), .combine=c) %dopar% {
    # how many p-value are lower among the random matrix (lower is better)
    length(which(p.rand <= p[i])) / R
}
qval <- qval/length(qval)

## COMBINING THE RESULTS TOGETHER
library('hgu133a2.db')
library(multtest)
e <- as.data.frame(e)
e <- cbind(pval=p, qval=qval, symbol=as.vector(unlist(mget(rownames(e), hgu133a2SYMBOL, ifnotfound=NA))), e)
e <- e[order(e$pval),]
e <- cbind(adj.pval=mt.rawp2adjp(e$pval)$adjp[,"BH"], e)

# How many gene have a permutation q-value <= 0.05?
dim(subset(e, qval<=0.05))
[1] 184 21
# How many gene have a Benjamini-Hochberg q-value <= 0.05?
dim(subset(e, adj.pval<=0.05))
[1]  0 22