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



t.differential <- function(
	x,					#signal data
	xcols,					#col numbers group 1
	ycols,					#col numbers group 2
	gin,					#genome information file , including "Gene.Symbol","label"
	label.col="label",		#name of gene description col
	parametric=T,			#filter by t test else by U test
	Shapiro.W.alpha=0.05,	#decision if normal distribution : >alpha: H0: normal distribution -> set flag, !:outlier sensitive
	F.p.alpha=0.05,			#decision if t-test or Welch variant : F test alpha error, >alpha: H0: no difference in var
	t.p.alpha=0.05,			#filter: t test alpha error, <=alpha: Ha: difference in means
	U.p.alpha=0.05,			#alternative for t test: non parametric: Wilcoxon Mann Whitney U test alpha error, <=alpha: Ha: difference in means
	F.filter=F,				#filter also according to var: keep only: >alpha: H0: no difference in var
	fold.change.k=1.5,		#filter: ratio
	diff.k=10,				#filter: difference of group means
	sig.file=F,				#save results in file
	sig.df=F,				#save results in data frame
	data.name="output",		#(name for output data frame(s))  &  graphics
	relative.path="/results/",	#path relative according to getwd() for graph sheets + text file
	get.all=F,				#get all calculated values - debug algorithm
	sort="no",				#sort output by "t.p.value", "difference"
	test.corr="BH",			#multiple testing correction  "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
	plot.it=T,				#plot graphs? - always be printed in a file -
	volcano=F,				#plot volcano graph?
	MARGIN=2,				#margins of plot
	rows.on.page=28,		#adjust number of result rows per plot page
	cex=0.6,				#text size
	lwd=3,					#line width
	font=10,					#font number
	color.set="greenBlackRed",	#change the color palette - see color.set()
	swscr=T					#T: standard reporting in file , F: to console (or triggered by other functions)
	)
{
	# compare microarray data, two groups
	# compute mean, t-test, confidence limits, differences, fold changes, F test values
	# Version 07022011
	
	# source dependencies
	source("/home/korschi/eclipseR/0functions/tdifferential/0.source.R")
	
	# checks 1
	if(missing(gin)){ cat("\n GIN missing \n"); return() }
	if(missing(xcols) | missing(ycols)){ cat("\n Group selection missing \n"); return() }
	if(length(xcols) < 2 | length(ycols) < 2){ cat("\n Size of each group must be > 1 \n"); return() }

	# ini
	require(genefilter)	#rowVars
	#Attaching package: The following objects are masked from ‘package:Biobase’:   anyMissing, rowMedians

	color.set(color.set)

	xname <- deparse(substitute(x))
	cxname <- names(x[,xcols])
	cyname <- names(x[,ycols])
	ncx <- length(xcols)
	ncy <- length(ycols)
	nr1 <- nrow(x)
	nc1 <- ncol(x)
 	rnames <- dimnames(x)[[1]]
	nxcols <- 1:ncol(x)

	# build descriptive statistics
	xcols.var <- rowVars(x[, xcols])		# might be 0
	xcols.rowMeans <- rowMeans(x[, xcols])
	ycols.var <- rowVars(x[, ycols])
	ycols.rowMeans <- rowMeans(x[, ycols])
	diff.means <- xcols.rowMeans - ycols.rowMeans
	fold.change <- xcols.rowMeans/ycols.rowMeans
	
	# handling of divison by 0  --  intercept Inf 1/0 or NaN 0/0  --  replaced by arbitary values
#	cat("\n inf ",sum(is.infinite(fold.change))," nan ",sum(is.nan(fold.change)))
	fold.change[is.infinite(fold.change)] <- 10
	fold.change[is.nan(fold.change)] <- 1
	# give warning
	cat("\n warning: fc values Inf/NaN detected and replaced by arbitary values 10 rsp. 1")
	
	# transform values <1 in: -1/x
	fold.change[fold.change<1] <- -1/fold.change[fold.change<1]
	
    # Shapiro-Wilk test (data normal distributed? - sample size 3-5000 !)
	if(ncx>2 & ncy>2){
		SW.result <- shapiro.wilk.test(x[,c(cxname,cyname)],
				#normal distribution = T , true for g1 and g2, AND ratio: W1/W2 deviation between the distributions
				Shapiro.W.alpha=Shapiro.W.alpha,
				ncx=ncx,
				ncy=ncy)
	}else{
		SW.result <- data.frame(matrix(NA,nr1,3))
		names(SW.result) <- c("sw.n1", "sw.n2", "sw.w.ratio")
	}
	
    # F-test
    F.result <- F.test.fast(rownames.x=rnames,
		#p value AND logical vector var equal = T, non equal = F
		ncx=ncx,
		ncy=ncy,
		xcols.var=xcols.var,
		ycols.var=ycols.var,
		alpha=F.p.alpha,
		alternative="two.sided")

    # t-test (two sided, two sample, with Welch - based on F-test results <=alpha: Ha: Var not equal, >alpha: H0: Var equal)
    t.result <- t.test.fast(rownames.x=rnames,
		var.equal=F.result[,2],               #var equal = T , true
		ncx=ncx,
		ncy=ncy,
		xcols.var=xcols.var,
		ycols.var=ycols.var,
		diff.means=diff.means,
		mu=0,
		alpha=t.p.alpha)

    # Wilcoxon rank sum test OR Mann-Whitney U test OR U test (non parametric test for normal distribution = F)
    U.result <- Wilcoxon.Mann.Whitney.U.test(x[,c(cxname,cyname)],
        #Wilcocon rank sum test: mu=0: H0: >alpha: mu1=mu2, Ha: <=alpha: mu1 unequal mu2
		ncx=ncx,
		ncy=ncy)

	rawtp <- t.result[,1]
	rawUp <- U.result[,1]
	
	# correction for multiple testing FWER/FDR			(private version: adjust.p)
	#  can take NA and NaN values
	if(test.corr!="none"){
		t.p.value <- p.adjust(
						p=rawtp,
						method=test.corr,
						n=length(rawtp))
		U.p.value <- p.adjust(
						p=rawUp,
						method=test.corr,
						n=length(rawUp))
	}else{
	    t.p.value <- rep(0,nr1)
	    U.p.value <- rep(0,nr1)
	}
	
    # build data structure of all results
	if(test.corr!="none"){
	    SIG <- factor(cut(t.p.value, breaks = c(0, 0.0001, 0.001, 0.01, 0.05, 100), labels = c("****", "***", "**", "*", "n"), include.lowest = T))
    }else{
	    SIG <- factor(cut(rawtp, breaks = c(0, 0.0001, 0.001, 0.01, 0.05, 100), labels = c("****", "***", "**", "*", "n"), include.lowest = T))
    }
	
	spec <- data.frame(cbind(xbar = xcols.rowMeans, ybar = ycols.rowMeans), row.names = rnames, check.rows = F)
	Gene.Symbol <- as.character(gin[rnames,"Gene.Symbol"])		#id and Gene.Symbol in result set
	diff.means <- as.double(diff.means)
	fold.change <- as.double(fold.change)
	rawtp <- as.double(rawtp)
	rawUp <- as.double(rawUp)
	ci <- as.double(t.result[,3])
	t.p <- as.double(t.p.value)     #corrected p
	U.p <- as.double(U.p.value)      #corrected p
	SIG <- as.character(SIG)
	sw.n1 <- SW.result[,1]
	sw.n2 <- SW.result[,2]
	sw.wratio <- SW.result[,3]
	rawFp <- F.result[,1]
	spec <- cbind(spec,Gene.Symbol,diff.means,fold.change,rawtp,rawUp,ci,t.p,U.p,SIG,sw.n1,sw.n2,sw.wratio,rawFp, stringsAsFactors=F)
	
	# open file for reporting
	if(swscr){ report(name=paste(relative.path,data.name,".report.",format(Sys.time(), "%Y%m%d%H%M"),".txt",sep="")) }    #switch from standard output to report file
	
	cat("\nSignificance Analysis of Microarray Data using t-Tests")
	cat("\nRows in x :",nr1)
	cat("\nData:",xname)
	cat("\nSummary of X:\n")
	print(summary(x))
	cat("\n")
	
	main <- paste("t-test comparisons with ", round((1-t.p.alpha)*100), "% confidence limits")
	if(test.corr!="none"){ main <- paste(main,", p-adjust.:",test.corr) }
	main <- paste(main,", data: ",xname,sep="")
	main <- paste(main,"\n Group 1:",sep="")
	for(i in 1:ncx){ main <- paste(main,cxname[i],sep=" ") }
	main <- paste(main,"\n Group 2:",sep="")
	for(i in 1:ncy){ main <- paste(main,cyname[i],sep=" ") }
	
	grp <- paste("\nSignificance Analysis of Microarray Data using t-Tests")
	grp <- paste(grp,"\nData:\t\t",xname)
	grp <- paste(grp,"\nRows:\t\t",nr1)
	grp <- paste(grp,"\nCols:\t\t",nc1)
	grp <- paste(grp,"\nt-test:")
	grp <- paste(grp,"\nGroup 1:\t",sep="")
	for(i in 1:ncx){ grp <- paste(grp,cxname[i],sep=" ") }
	grp <- paste(grp,"\nGroup 2:\t",sep="")
	for(i in 1:ncy){ grp <- paste(grp,cyname[i],sep=" ") }
	grp <- paste(grp,"\nt alpha:\t",t.p.alpha)
	if(test.corr!="none")grp <- paste(grp,"\np-adjustment:\t",test.corr)
	
	cat("\n",main)
	cat("\n",grp)
	
    # build HO subset step by step
	sel.set <- rep(T,nr1)  #running sum
	sel.filter <- sel.set  #indicator for filter performance
	sel.mem <- sel.set     #regenerate filter vector
	
	cat("\n\nComparisons and cutsets for different tests and threshold levels")
	cat("\n  N      Test            Threshold         Intersection")
	cat("\n________________________________________________________________")
	cat("\n",sum(sel.set),"\tComplete H0")
	
    if(parametric){
         sel.set[rawtp>t.p.alpha] <- F
         sel.filter[rawtp>t.p.alpha] <- F
         cat("\n",sum(sel.filter),"\tt-test            p level <=",t.p.alpha,"\t",sum(sel.set))
         sel.filter <- sel.mem

	     if(test.corr!="none"){
             sel.set[t.p.value>t.p.alpha] <- F
             sel.filter[t.p.value>t.p.alpha] <- F
             cat("\n",sum(sel.filter),"\tadjusted t p-values, method:",test.corr,"\t",sum(sel.set))
		     sel.filter <- sel.mem
         }
    }else{
         sel.set[rawUp>U.p.alpha] <- F
         sel.filter[rawUp>U.p.alpha] <- F
         cat("\n",sum(sel.filter),"\tU-test            p level <=",U.p.alpha,"\t",sum(sel.set))
         sel.filter <- sel.mem

	     if(test.corr!="none"){
             sel.set[U.p.value>U.p.alpha] <- F
             sel.filter[U.p.value>U.p.alpha] <- F
             cat("\n",sum(sel.filter),"\tadjusted U p-values, method:",test.corr,"\t",sum(sel.set))
		     sel.filter <- sel.mem
         }
    }
	
    if(F.filter){
        sel.set[rawFp<F.p.alpha] <- F            #keep only those with Var equal in both groups
        sel.filter[rawFp<F.p.alpha] <- F
        cat("\n",sum(sel.filter),"\tF-test            p level <=",F.p.alpha,"\t",sum(sel.set))
        sel.filter <- sel.mem
    }
	
    sel.set[(fold.change<fold.change.k) & (fold.change>(-fold.change.k))] <- F
    sel.filter[(fold.change<fold.change.k) & (fold.change>(-fold.change.k))] <- F
	cat("\n",sum(sel.filter),"\tFold Change           fc  >=",fold.change.k,"\t",sum(sel.set))
    sel.filter <- sel.mem
    
    sel.set[(diff.means<diff.k) & (diff.means>-diff.k)] <- F
    sel.filter[(diff.means<diff.k) & (diff.means>-diff.k)] <- F
	cat("\n",sum(sel.filter),"\tDifference of group means   >=",diff.k,"\t",sum(sel.set))
    sel.filter <- sel.mem

    cat("\n\n\tFinal intersection              \t",sum(sel.set),"\n")

	cat("\n________________________________________________________________\n\n")

    # sort it - if desired
	if(sort=="t.p.value"){
		op <- order(spec$t.p.value)
		spec <- spec[op,]
	}
	if(sort=="difference"){
		xbar <- spec[, "xbar"]
		ybar <- spec[, "ybar"]
		od <- rev(order(abs(xbar-ybar)))
		spec <- spec[od,]
	}

	# output all in data.frame, plot all and exit
	# -- also good for looking at a small set of factors independently of their significance --
	if(get.all){
		fn.var <- paste(getwd(),relative.path,data.name,format(Sys.time(), "%Y%m%d%H%M"),".txt", sep="")
		cat("\n results in ",fn.var," \n")
		cat("\n ALL mode - All calculations saved \n")
		if(swscr){ script() }    #switch from report window to standard output
		
		if(plot.it){		#plot?
			plot.differential(
					x=spec,		#all values
					t.or.U=if(parametric){ "t" }else{ "U" },
					alpha=if(parametric){ t.p.alpha }else{ U.p.alpha },
					fold.change.k=fold.change.k,
					diff.k=diff.k,
					gin=gin,
					test.corr=test.corr,
					label.col=label.col,		
					main=main,
					rows.on.page=rows.on.page,
					cex=cex,
					lwd=lwd,
					data.name=data.name,
					relative.path=relative.path
			)
			
			if(volcano){
				plot.volc(
						x=spec,			#all values
						fc="fold.change",
						tp=if(test.corr=="none"){ if(parametric){ "rawtp" }else{ "rawUp" } }else{ if(parametric){ "t.p" }else{ "U.p" } },
						dm="diff.means",
						treshold.fc=fold.change.k,
						treshold.tp=if(parametric){ t.p.alpha }else{ U.p.alpha },
						cex=cex,
						data.name=data.name,
						relative.path=relative.path
				)
			}
		}
		
		# save all data
		if(sig.file){
			write.table(x=spec, file=fn.var,
					sep="\t", row.names=T, col.names=T, quote=F)
			cat("\n saved in ",fn.var,"\n")
		}else{
			cat("\n all.file=F \n")
		}
		if(sig.df){
			assign(x=data.name, value=spec, pos=1)
		}else{
			cat("\n all.df=F \n")
		}
		
		return()	#stop
    }
 	
	
	if(sum(sel.set)==0){	# no graphics because of no results
		cat("\n no results \n")
		if(swscr){ script() }    #switch from report file to standard output
		
		plot.volc(		# take a look on the data
				x=spec,			#all values
				fc="fold.change",
				tp=if(test.corr=="none"){ if(parametric){ "rawtp" }else{ "rawUp" } }else{ if(parametric){ "t.p" }else{ "U.p" } },
				dm="diff.means",
				treshold.fc=fold.change.k,
				treshold.tp=if(parametric){ t.p.alpha }else{ U.p.alpha },
				cex=cex,
				data.name=paste(data.name,".no.sig",sep=""),
				relative.path=relative.path
		)
		
		return(0)	#stop processing
	}
	
    Significant <- spec[sel.set,,drop=F]	# create result set
	
	fn.var <- paste(getwd(),relative.path,data.name,".",format(Sys.time(), "%Y%m%d%H%M"),".txt", sep="")
	cat("\n results in ",fn.var," \n")
	
	if(swscr){ script() }    #switch from report file to standard output

    # plot results
	if(plot.it){		#commented plot 1?
		plot.differential(
				x=Significant,		#only significant values
				t.or.U=if(parametric){ "t" }else{ "U" },
				alpha=if(parametric){ t.p.alpha }else{ U.p.alpha },
				fold.change.k=fold.change.k,
				diff.k=diff.k,
				gin=gin,
				test.corr=test.corr,
				label.col=label.col,		
				main=main,
				rows.on.page=rows.on.page,
				cex=cex,
				lwd=lwd,
				data.name=data.name,
				relative.path=relative.path
		)

		if(volcano){
			plot.volc(
				x=spec,			#all values
				fc="fold.change",
				tp=if(test.corr=="none"){ if(parametric){ "rawtp" }else{ "rawUp" } }else{ if(parametric){ "t.p" }else{ "U.p" } },
				dm="diff.means",
				treshold.fc=fold.change.k,
				treshold.tp=if(parametric){ t.p.alpha }else{ U.p.alpha },
				cex=cex,
				data.name=data.name,
				relative.path=relative.path
			)
		}
	}
	
	#
	palette("default")
	
	if(sig.file){
		write.table(x=Significant, file=fn.var, sep="\t", row.names=T, col.names=T, quote=F)
		cat("\n saved in ",fn.var,"\n")
	}else{
		cat("\n sig.file=F \n")
	}
	if(sig.df){
		assign(x=data.name, value=Significant, pos=1)
	}else{
		cat("\n sig.df=F \n")
	}
	return()
}