Skip to content

Instantly share code, notes, and snippets.

@bobthecat
Last active December 9, 2015 23:58
Show Gist options
  • Save bobthecat/4347526 to your computer and use it in GitHub Desktop.
Save bobthecat/4347526 to your computer and use it in GitHub Desktop.
## 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment