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



# use house keeping genes / ERCC spike in controls


normHouseKeeping <- function(x, y, hk=c("ACTB","GAPDH","PGK1"), ref.col=1, ps="mean", subfolder.filename=NULL){
	# normalize according to some houskeeping genes which are part of the input data
	# x: data to be normalized - columns: experiments, rows: genes / exons ... (data.frame)
	# y: gene symbol column (corresponding row order to x, vector)
	# hk: list with house keeping genes: gene symbols - somthing which can be found
	#  in the y data vector
	# ref.col: column which remains to be unchanged
	# ps: if more than one reference gene row: use mean or median to produce an column
	#  average to be used as a normalization factor
	# subfolder.filename: path + name for results
	
	#ini
	nrx <- nrow(x)
	ncx <- ncol(x)
	nry <- length(y)
	if(nrx!=nry){ cat("\n length of x and y are not allowed to differ"); break }
	if(ncx<2){ cat("\n minimum of 2 columns in x"); break }
	len.hk <- length(hk)
	hits <- NULL
	
	#work
	#find lines
	for(i in 1:len.hk){
		hits <- c( hits, grep(pattern=paste("^",hk[i],"$",sep=""), x=y,
				ignore.case = FALSE, perl = FALSE, value = FALSE, fixed = FALSE, useBytes = FALSE, invert = FALSE) )
	}
	
	#select subset
	sel.lines <- x[hits,]
	#order ref col on first position and all other thereafter
	sel.lines.o <- sel.lines[, c(ref.col, c(1:ncx)[-ref.col]) ]
	sel.lines.o.nr <- nrow(sel.lines.o)		#count rows
	#caculate transformation factors for each row
	f1 <- function(x){
		len.x <- length(x)
		erg.factor <- vector(mode="numeric",length=len.x-1)
		for(i in 2:len.x){
			erg.factor[i-1] <- x[1] / x[i]
		}
		return(erg.factor)
	}
	erg.factor <- t( apply(sel.lines.o,1,f1) )
	#check for NaN (0/0) or NA (0/NA)
	erg.factor[is.na(erg.factor)] <- 0
	#check for Inf (1/0)
	erg.factor[is.infinite(erg.factor)] <- 0
	#return(erg.factor)
	
	#plot distribution per factor
	for(i in 1:(ncx-1)){
		hist.plot.simple(x=erg.factor[,i],
			bin.num=F,bin.size=F,xlab="",xlim=NULL,y.max=NULL,y.at=NULL,y.num.format=F,ylog=F,
			bar.width=1,offset=0,digits=1,col=c("red","green","blue"),cex=1,norm=F,add=F,n.factor=1,curve=F,smooth.f=0,
			subfolder.filename=paste(subfolder.filename, ".", "HK", i, ".", ps, ".", deparse(substitute(x)), sep="") )
	}
	
	#construct normalization factor per column
	if(sel.lines.o.nr>1){
		if(ps=="mean"){
			erg.factor.m <- apply(erg.factor, 2, mean, na.rm=FALSE)
		}
		if(ps=="median"){
			erg.factor.m <- apply(erg.factor, 2, median, na.rm=FALSE)
		}
	}else{
		erg.factor.m <- erg.factor
	}
	
	#show values
	if(ref.col==1){ cat("\n Factor(s) ref ", erg.factor.m, "\n") }
	if(ref.col>1 & ref.col<ncx){
		seq1 <- (1:(ref.col-1))
		seq2 <- (ref.col:(ncx-1))
		cat("\n Factor(s) ", erg.factor.m[seq1], " ref ", erg.factor.m[seq2], "\n")
	}
	if(ref.col==ncx){ cat("\n Factor(s) ", erg.factor.m, " ref \n") }
	
	
	#apply normalization factor per column
	seq1 <- (1:ncx)[-ref.col]
	j <- 0
	for(i in seq1){
		j <- j+1
		x[,i] <- x[,i]*erg.factor.m[j]
	}
	
	#
	return(round(x,digits=2))
}



normSpikeIn <- function(x, y, ref.col=1, factor.type="mean", subfolder.filename=NULL){
	# normalize according to ERCC spike in controls
	# x: data to be normalized - cols: experiments, rows: genes / exons ... (data.frame)
	# y: normalization data set (corresponding col order to x!, data.frame)
	# ref.col: column which remains to be unchanged
	# factor.type: if more than one row in y: use mean or median to produce an column
	#  average to be used as a normalization factor
	# subfolder.filename: path + name for results
	
	#ini
	nrx <- nrow(x)
	ncx <- ncol(x)
	nry <- nrow(y)
	if(ncx<2){ cat("\n minimum of 2 columns in x"); break }
	
	#work
	#order ref col on first position and all other thereafter
	y.o <- y[, c(ref.col, c(1:ncx)[-ref.col]) ]
	
	#caculate transformation factors for each row
	f1 <- function(x){
		len.x <- length(x)
		erg.factor <- vector(mode="numeric",length=len.x-1)
		for(i in 2:len.x){
			erg.factor[i-1] <- x[1] / x[i]
		}
		return(erg.factor)
	}
	erg.factor <- t( apply(y.o,1,f1) )
	#check for NaN (0/0) or NA (0/NA)
	erg.factor[is.na(erg.factor)] <- 0
	#check for Inf (1/0)
	erg.factor[is.infinite(erg.factor)] <- 0
	assign("detailsnormSpikeIn", value=erg.factor)		#write in present workspace
	
	#plot distribution per factor
	for(i in 1:(ncx-1)){
		hist.plot.simple(x=erg.factor[,i],
				bin.num=F,bin.size=F,xlab="",xlim=NULL,y.max=NULL,y.at=NULL,y.num.format=F,ylog=F,
				bar.width=1,offset=0,digits=1,col=c("red","green","blue"),cex=1,norm=F,add=F,n.factor=1,curve=F,smooth.f=0,
				subfolder.filename=paste(subfolder.filename, ".", "Spike", i, ".", factor.type, ".", deparse(substitute(x)), sep="") )
	}
	
	#construct normalization factor per column
	if(nry>1){
		if(factor.type=="mean"){
			erg.factor.m <- apply(erg.factor, 2, mean, na.rm=FALSE)
		}
		if(factor.type=="median"){
			erg.factor.m <- apply(erg.factor, 2, median, na.rm=FALSE)
		}
	}else{
		erg.factor.m <- erg.factor
	}
	
	#show values
	if(ref.col==1){ cat("\n Factor(s) ref ", erg.factor.m, "\n") }
	if(ref.col>1 & ref.col<ncx){
		seq1 <- (1:(ref.col-1))
		seq2 <- (ref.col:(ncx-1))
		cat("\n Factor(s) ", erg.factor.m[seq1], " ref ", erg.factor.m[seq2], "\n")
	}
	if(ref.col==ncx){ cat("\n Factor(s) ", erg.factor.m, " ref \n") }
	
	
	#apply normalization factor per column
	seq1 <- (1:ncx)[-ref.col]
	j <- 0
	for(i in seq1){
		j <- j+1
		x[,i] <- x[,i]*erg.factor.m[j]
	}
	
	#
	return( round(x,digits=2) )
}