# TODO: Add comment
# 
# Author: E.Korsching 10.9.2009
###############################################################################



plot.volc <- function(
	x,						#data.frame containing the following cols
	fc,						#parameter 1: col name : e.g. "fold.change" - input range: -n..-1,1..n
	tp,						#parameter 2: col name : e.g. "t.p.value"
	dm=NULL,				#parameter 3: col name : e.g. "diff.means"		(parameter optional)
	treshold.fc=NULL,		#treshold zu 1 : values : e.g. fold.change.k	(parameter optional)
	treshold.tp=NULL,		#treshold zu 2 : values : e.g. t.p.alpha		(parameter optional)
	cex=1,					#text size
	data.name="test",		#plot in file
	relative.path="/results/"
)
{
	# Volcano plot : expression ratio (fold change) log2(mean 1 / mean 2)  vs. statistical significance e.g. t-tests [-log10(p)] 
	# ref.: Welle S, Brooks AI et al., Computational method for reducing variance with Affymetrix microarrays, BMC Bioinformatics 2002,3:23
	# cols c1, c2, c3 have to be the same length
	# highlight: e.g. points(log2(ratio[interesting]), -log10(p[interesting]), pch=16, col="red")
	
	# ini
	#colors: set 1..16 ramp before e.g. by color.set()
	nrx <- nrow(x)
	xname <- deparse(substitute(x))
	rnames <- row.names(x)
	
	#transform fold change values back for 'volcano' purpose : -n..-1 with -1/x,  1..n values remain unchanged
	fc.logical <- x[,fc]<0
	fc.x <- x[,fc]
	fc.x[fc.logical] <- -1/x[fc.logical,fc]
	
	#transform x,y values according to log2 and log10 scale
	fc.log2.x <- log2(fc.x)
	treshold.fc.log2.x <- log2(treshold.fc)
	tp.log10.y <- -log10(x[,tp])
	treshold.tp.log10.y <- -log10(treshold.tp)
	
	#	if(test.corr!="noadjust"){ ylimit<- c(0,abs(logb(t.p.value/100,10))) }else{ ylimit<- c(0,abs(logb(rawtp/100,10))) }
	
	# open graph device
#	png(filename=paste(getwd(),relative.path,data.name,".volc.",format(Sys.time(), "%Y%m%d%H%M."),".png", sep=""),
#			width=800, height=800, units="px", pointsize=12, bg="white")
	pdf(file=paste(getwd(),relative.path,data.name,".volc.",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""),
			width=11, height=7, onefile=T, title=paste(data.name,".volcano plot", sep=""), pointsize=12)
	
	par(mar=c(8, 5, 5, 5), cex=cex)
	plot(fc.log2.x,  tp.log10.y, xlab="", ylab="", pch=1, axes=F)		# xlim=xlimit, ylim=ylimit,
	xy.usr <- par("usr")			# get the created usr coordinates
	
	#colorize points by c3 range
#	if(!is.null(dm)) {
#		dm.linear <- x[,dm]
#		#translate dm values to color values according to palette(16)
#		dm.range <- range(dm.linear)
#		dm.steps <- seq(dm.range[1], dm.range[2], (dm.range[2]-dm.range[1])/16)
#		z <- dm.linear
#		for(i in 1:16){
#			z[dm.linear<=dm.steps[i]] <- i
#		}
#		#
#		for(i in 1:16) {			# loop through all colors
#			points(fc.log2.x[z==i], tp.log10.y[z==i], pch=16, col=palette()[i], cex=3*cex)
#	}
	
	#or highlight values by thresholds
#	if(is.null(dm) & !is.null(treshold.fc) & !is.null(treshold.tp) ){
		# -treshold.fc.log2.x  < x <  treshold.fc.log2.x
		tmp.logic.a <- fc.log2.x>treshold.fc.log2.x | fc.log2.x< -treshold.fc.log2.x
		#treshold.tp.log10.y
		tmp.logic.b <- tp.log10.y>treshold.tp.log10.y	# > : because of: -log approach
		#
		tmp.logic <- tmp.logic.a & tmp.logic.b
		#cat("\n tmp.logic.a ",tmp.logic.a)
		#cat("\n tmp.logic.b ",tmp.logic.b)
		#cat("\n tmp.logic   ",tmp.logic)
		points(fc.log2.x[ tmp.logic ], tp.log10.y[ tmp.logic ], pch=16, col="green", cex=3*cex)
#	}
	
	#creating values for axis, grid
	xd <- signif(xy.usr[2]-xy.usr[1], 1)
	xs <- signif(xy.usr[1]/xd, 1)
	xa <- xd*xs
	xd.t <- signif(xd/10, 1)
	xx <- seq(xa,xy.usr[2],xd.t)	#x custom tick marks
	yd <- signif(xy.usr[4]-xy.usr[3], 1)
	ys <- signif(xy.usr[3]/yd, 1)
	ya <- yd*ys
	yd.t <- signif(yd/10, 1)
	yy <- seq(ya,xy.usr[4],yd.t)	#y custom tick marks
	
	yy3 <- 2^xx			# show true (de-log2) x values
	yy4 <- 10^(-yy)	# show true (de-log10) y values
	axis(1, at=xx)
	axis(2, at=yy)
#	axis(1)
#	axis(2)
	axis(3, at=xx, labels=round(yy3, 2))
	axis(4, at=yy, labels=round(yy4, 2), adj=0)
#	abline(v=xx, lty=1, lwd=0.2, col="black")
#	abline(h=yy, lty=1, lwd=0.2, col="black")
	abline(v=c(treshold.fc.log2.x, -treshold.fc.log2.x), lty=2, lwd=0.3, col="black")		# +/- fc limit
	abline(h=treshold.tp.log10.y, lty=2, lwd=0.3, col="black")				# significance limit
	
	## axis legend in the margins
	for(i in 1:4){
		mtext(text=c("log2(fold change)","-log10(p value)","fold change","p value")[i], side=i, line=3, font=2)
	}
	
	#draw box around plot
	box("figure")
	
	#figure text
	sub <- paste("\nVolcano plot : \'", xname, "\' data set", sep="")
	title(sub=sub, adj=0)
	
#	if(!is.null(dm)) {			#optional legend to plot parameter c3
#		dm.linear <- x[,dm]	#only test -> del after
#		lab <- round(seq(min(dm.linear), max(dm.linear), length=16))
#		legend(x=xy.usr[1]+(0.1*(xy.usr[2]-xy.usr[1])), y=(xy.usr[2]-xy.usr[1])/2, legend=lab, density=-1, fill=palette(), yjust=1, border="black", bty="n", cex=cex)
#	}
	
	# close graph device
	dev.off()
}