analyse_layerEnrich.R

AJWillsey — Mar 25, 2013, 2:24 PM

rm(list=ls())
setwd("~/Dropbox/SSC_July_2012/Co-expression/coexp_paper_Nov-12/newest_analysis/03-25-2013_LayerEnrich/")

#read in file
mydata = read.table("all_dex.txt", head=F, sep="\t");
#mydata=mydata[,-1]

#define function for getting cumulative sum of each number of hits; -log10(p-value) of number of hits
myCumFunctionSimple=function(x){
  result=rev(cumsum(rev(table(x)))) 
  return(result)
}

myLogFunction=function(x){
  result=(-1)*log10(x/20000)
  return(result)
}

#apply functions
cumdata=as.data.frame(apply(mydata, 2, myCumFunctionSimple), ncol=1)
rownames(cumdata)=rownames(apply(mydata, 2, myCumFunctionSimple))
logData=apply(cumdata,2, myLogFunction)
rownames(logData)=rownames(cumdata)

#plot curve
#setEPS()
#postscript(file="all_all_permutations.eps", width=8.5, height=3, colormodel='cmyk')
#par(mfcol=c(1,2), mai=c(0.75, 1, 0.25 ,0.5), omi=c(0,0,0,0), ps=12, font.lab=2, font.main=2,
#    mgp=c(2, 0.5, 0), las=1, cex=1, cex.lab=1, cex.main=1, cex.axis=1)

plot(lowess(as.numeric(rownames(logData)), logData), ylim=c(0,5), xlab="CPi + CPo + SP enriched genes", ylab="-log10(p)", main="Significance of layer-enriched genes", pch=21, xaxs='r', type='l', lwd=2, col="black")
abline(h=1.301, col="#922686", lty=2, lwd=2)
abline(v=14, col="#ce2a37", lwd=2)

legend("topleft", inset=0.025, c("Permutation", "Observed", "p=0.05"),
       lty=c(1, 1, 2), lwd=c(2, 2, 2), col=c("black", "#ce2a37", "#922686"))
#legend("topleft", inset=0.025, c("Permutation", "Observed"), lty=c(0,1), lwd=c(0,2), pch=c(22, NA), pt.bg=c("#6FCCDD", "#ce2a37"), pt.lwd=1.5, col=c("black", "#ce2a37"), pt.cex=c(2,2)) 

mtext(paste("p=", signif(cumdata[which(rownames(cumdata)=="14"),]/length(mydata[,1]), digits=2), sep=""), side=3)

#dev.off()