R/MultiRankSeqFunctions.R

In slzhao/MultiRankSeq: MultiRankSeq

# TODO: Add comment
# 
# Author: zhaos
###############################################################################

#' @import DESeq2 edgeR baySeq VennDiagram knitr
#' @name exampleCounts
#' @title An example data generated by MultiRankSeq package
#' @description This data set is a random counts data based on negative binomial distribution. The codes for generating this data can be found on the example of findDiff function.
#' @docType data
#' @usage exampleCounts
NULL
#' @name exampleCountsDiff
#' @title An example data generated by MultiRankSeq package
#' @description This data set is the results from DESeq, edgeR, baySeq packages, which is based on the exampleCounts data set. The codes for generating this data can be found on the example of findDiff function.
#' @docType data
#' @usage exampleCountsDiff
NULL
##' MultiRankSeqReport
##' 
##' A wrapper for function findDiff and diffReport.
##' 
##' A function to generate MultiRankSeq report based on counts data.
##' 
##' @param output the path and name for report file.
##' @param rawCounts the count data used in analysis, which should be a matrix or a data.frame.
##' @param group the group of count data, which should be 0 or 1. The samples labeled as 1 will be compared with samples labeled as 0.
##' @param paired the pairs of count data. Will be NULL if no pair. For example, if sample 1 and 2 are one pair, sample 3 and 4 are another pair, the paired should be c(1,1,2,2).
##' @param seed the random seed for bayseq. The user can set it so that the bayseq result can be reproduced.
##' @param diffResult The result generated by \code{\link{findDiff}} function. If not specified, findDiff function will be used to generate it.
##' @param exportCounts Logical. Indicate if the raw counts will be exported in result table.
##' @param useVStCount Logical. The DESeq2 package can apply a variance stabilizing transformation (VST) to count data, which can be used in checking outliers or clustering. This parameter indicates if the transformed counts will be used in heatmap, boxplot, and correlation in reports.
##' @param useAllInHeatmap Logical. 
##' @param percent The percent of genes with highest SD, to be used in heatmap. 
##' @param ... Other parameters used in \code{\link{diffReport}} function.
##' @export
##' @examples #A simple example
##' data(exampleCounts)
##' data(exampleCountsDiff)
##' result<-MultiRankSeqReport(output="MultiRankSeq1.html",diffResult=exampleCountsDiff,rawCounts=exampleCounts)
##' \dontrun{
##' #An example: from counts data to report, paired data
##' #this code may take two minutes, because performing baySeq, DESeq2, and edgeR may be slow. 
##' result<-MultiRankSeqReport(output="MultiRankSeq2.html",rawCounts=exampleCounts,group=c(0,0,0,1,1,1),paired=c(1:3,1:3))
##' }
MultiRankSeqReport<-function(output="MultiRankSeq.html",rawCounts,group=c(0,0,0,1,1,1),diffResult,seed=NULL,paired=NULL,exportCounts=F,useVStCount=T,useAllInHeatmap=F,percent=0.05,...) {
	rawCounts<-rawCounts[which(rowSums(rawCounts)>0),]
	if (missing(diffResult)) {
		diffResult<-findDiff(rawCounts=rawCounts,group=group,paired=paired,seed=seed)
	}
	diffReport(output=output,diffResult=diffResult,rawCounts=rawCounts,group=group,paired=paired,exportCounts=exportCounts,useVStCount=useVStCount,useAllInHeatmap=useAllInHeatmap,percent=percent,...)
	return(diffResult)
}

##' findDiff
##' 
##' The function will perform DESeq, edgeR, and baySeq analysis.
##' 
##'
##' 
##' @inheritParams MultiRankSeqReport
##' @export
##' @return A matrix with log2 fold change, p value, adjusted p value, and rank of DESeq, edgeR, and baySeq analysis.
##' @examples set.seed(123456789) #generate data
##' ngenes <- 500 #1000
##' q0 <- rexp(ngenes, rate = 1/250)
##' is_DE <- runif(ngenes) < 0.3
##' lfc <- rnorm(ngenes, sd = 2)
##' q0A <- ifelse(is_DE, q0 * 2^(lfc/2), q0)
##' q0B <- ifelse(is_DE, q0 * 2^(-lfc/2), q0)
##' true_sf <- c(1, 1.3, 0.7, 0.9, 1.6,1.5)
##' conds <- c("A", "A", "A","B", "B", "B")
##' exampleCounts <- t(sapply(seq_len(ngenes), function(i) sapply(1:6, function(j) rnbinom(1,mu = true_sf[j] * ifelse(conds[j] == "A", q0A[i], q0B[i]),	size = 1/0.2))))
##' colnames(exampleCounts)<-c(paste("Control",1:3,sep=""),paste("Disease",1:3,sep=""))
##' row.names(exampleCounts)<-paste("Gene",1:ngenes,sep="")
##' #this code may take two minutes, because performing baySeq, DESeq2, and edgeR may be slow. 
##' exampleCountsDiff<-findDiff(exampleCounts) #find diff result
findDiff<-function(rawCounts,group=c(rep(0,ncol(rawCounts)/2),rep(1,ncol(rawCounts)/2)),seed=NULL,paired=NULL) {
	rawCounts<-rawCounts[which(rowSums(rawCounts)>0),]
	if (is.null(row.names(rawCounts))) {
	  row.names(rawCounts)=1:nrow(rawCounts)
	}
	if (is.null(colnames(rawCounts))) {
	  colnames(rawCounts)=1:ncol(rawCounts)
	}
	comparDESeq<-myDESeq2(rawCounts,group,paired)
	comparEdgeR<-myedgeR2(as.matrix(rawCounts),group,paired)
	if (!is.null(seed)) {
		set.seed(seed)
	}
	comparBayseq<-mybayseq2(as.matrix(rawCounts),group,paired)
	resultCombinedCompar<-combined_result2(DESeq=comparDESeq,edgeR=comparEdgeR,bayseq=comparBayseq,rawCount=rawCounts,comparGroups=group)
	return(resultCombinedCompar)
}

##' diffReport
##' 
##' The function will generate a report based on the result of DESeq, edgeR, and baySeq analysis.
##' 
##'
##' 
##' @inheritParams MultiRankSeqReport
##' @param input The templet file. You don't need to specify it in most cases. A file in MultiRankSeq package will be used.
##' @param diffResult The result generated by \code{\link{findDiff}} function.
##' @param percent The percent of genes with highest SD, to be used in heatmap. 
##' @export
##' @examples data(exampleCounts)
##' data(exampleCountsDiff)
##' diffReport(output="MultiRankSeq.html",diffResult=exampleCountsDiff,rawCounts=exampleCounts)
diffReport<-function(input=system.file("extdata", "report.Rhtml", package = "MultiRankSeq"),output,diffResult,rawCounts,group=c(0,0,0,1,1,1),paired=NULL,exportCounts=F,useVStCount=T,useAllInHeatmap=F,percent=0.05) {
#	require(knitr)
#	require(VennDiagram)
	outputWd<-dirname(output)
	if (outputWd != ".") { #file name and path
		oldWd<-getwd()
		setwd(outputWd)
	} else { #only file name
		output<-paste(getwd(),output,sep="/")
	}
	if (exportCounts) {
		write.csv(cbind(diffResult,rawCounts),paste(output,".diff",".csv",sep=""))
	} else {
		write.csv(diffResult,paste(output,".diff",".csv",sep=""))
	}
	knit(input=input,output=output,quiet=T)
	print (paste("The MultiRankSeq report was generated in ",output,sep=""))
	if (outputWd != ".") {
		setwd(oldWd)
	}
}

plot3D<-function(diffResultPlot,...) {
	if (length(which(is.na(diffResultPlot[,1]) | is.infinite(diffResultPlot[,1])))>0) {
		diffResultPlot<-diffResultPlot[-which(is.na(diffResultPlot[,1]) | is.infinite(diffResultPlot[,1])),]
	}
	col<-rep("grey",nrow(diffResultPlot))
	col[which(diffResultPlot[,1]<=-1& diffResultPlot[,2]<=0.05)]<-"green"
	col[which(diffResultPlot[,1]>=1& diffResultPlot[,2]<=0.05)]<-"red"
	cex<-rep(1,nrow(diffResultPlot))
	temp<-1/100*(101-diffResultPlot[,3])
	temp[temp<0]<-0
	cex<-cex+temp*2
	plot(diffResultPlot[,1],-log10(diffResultPlot[,2]),pch=16,col=col,cex=cex,xlab="log2 fold change",ylab="-log10 p value",...)
	legend("top",legend=c("Rank=1","Rank>100"),pt.cex =c(3,1),pch=16,col="grey")
}

mybayseq2<-function(count, group,paired=NULL){
#	require(baySeq)
	cl <- NULL
	if (is.null(paired)) {
		conds <- list(NDE = rep(1, length(group)), DE = group)
		CD<-new("countData", data = count, replicates = group, groups = conds,
		        annotation=as.data.frame(row.names(count)))
#		CD@libsizes <- getLibsizes(CD, estimationType="quantile")
		libsizes(CD)<- getLibsizes(CD, estimationType="quantile")
		CDP.NBML <- getPriors.NB(CD, samplesize = 1000, 
				estimation = "QL", equalDispersions=TRUE, cl = cl)
#		CDPost.NBML <- getLikelihoods.NB(CDP.NBML, pET = 'BIC', cl = cl)
		CDPost.NBML <- getLikelihoods(CDP.NBML, pET = 'BIC', cl = cl)
		bayseq<-topCounts(CDPost.NBML, group = "DE", 
				number=dim(count)[1] )
		row.names(bayseq)=bayseq$row.names.count.
#		bayseq=bayseq[,-c(1:(dim(count)[2]+1))]
		if ("FDR.DE" %in% colnames(bayseq)) {
			bayseq<- bayseq[,c("likes","FDR.DE")]
		} else if ("FDR" %in% colnames(bayseq)) {
			bayseq<- bayseq[,c("likes","FDR")]
		} else {
			bayseq<- bayseq[,-c(1:(dim(count)[2]+1))][,c(1,ncol(bayseq))]
		}
	} else { #paired data
		if (length(which(group==0))*2!=length(group)) {
			stop("The two groups should have same number of samples for paired data") 
		}
		if (length(which(table(paired)!=2))) {
			stop("The paired data should have 2 samples for each pair") 
		}
		if (length(paired)!=length(group) | length(paired)!=ncol(count)) {
			stop("The paired, group, and raw count should have the same number of samples")
		}
		pairGroup<-sapply(unique(paired),function(x) which(paired==x))
		conds <- list(NDE = rep(1, (length(group)/2)))
		pairCD<-new("countData", 
				data = array(c(count[,pairGroup[1,]],count[,pairGroup[2,]]),dim = c(nrow(count), ncol(pairGroup), 2)),
				replicates = rep(1, (length(group)/2)), 
				groups = conds,
				densityFunction = bbDensity,
				annotation=as.data.frame(row.names(count)))
		libsizes(pairCD) <- getLibsizes(pairCD)
		pairCD <- getPriors(pairCD, samplesize = 1000, cl = cl)
		pairCD <- getLikelihoods(pairCD, pET = 'BIC', nullData = TRUE,cl = cl)
		bayseq<-topCounts(pairCD, group = 1,number=Inf)
#		row.names(bayseq)<-row.names(count)[bayseq$row.names.count.]
		# if ("FDR.NDE" %in% colnames(bayseq)) {
		# 	bayseq<- bayseq[,c("Likelihood","FDR.NDE")]
		# } else if ("FDR" %in% colnames(bayseq)) {
		# 	bayseq<- bayseq[,c("Likelihood","FDR")]
		# } else {
		# 	bayseq<- bayseq[,c("Likelihood",colnames(bayseq)[ncol(bayseq)])]
		# }
		row.names(bayseq)=bayseq$row.names.count.
		if ("FDR.NDE" %in% colnames(bayseq)) {
		  bayseq<- bayseq[,c("likes","FDR.NDE")]
		} else if ("FDR" %in% colnames(bayseq)) {
		  bayseq<- bayseq[,c("likes","FDR")]
		} else {
		  bayseq<- bayseq[,-c(1:(dim(pairGroup)[2]+1))][,c(1,ncol(bayseq))]
		}
	}
	if(colnames(bayseq)[1]=="likes") {
		bayseq[,1]<-1-bayseq[,1]
		colnames(bayseq)[1]<-"OneMinusLikelihood"
	}
	return(bayseq)
}

myDESeq2<-function(count, group,paired=NULL){
#	require(DESeq2)
	if (is.null(paired)) {
		dds <- DESeqDataSetFromMatrix(countData = count,
				colData = data.frame(group=as.factor(group)),
				design = ~group)
	} else {
		dds <- DESeqDataSetFromMatrix(countData = count,
				colData = data.frame(paired=as.factor(paired),group=as.factor(group)),
				design = ~paired+ group)
	}
	dds <- DESeq(dds)
	dds <- results(dds)
#	dds <- results(dds,independentFiltering=F)
	return(dds)
}

myedgeR2<-function(count, group,paired=NULL){
#	require(edgeR)
	y<-DGEList(count)
	group<-as.factor(group)
	if (is.null(paired)) {
		design <- model.matrix(~group)
	} else {
		paired<-as.factor(paired)
		names(paired)<-colnames(count)
		design <- model.matrix(~paired+group)
#		design <- model.matrix(~0+group+paired)
#		row.names(design)<-colnames(count)
	}
	y <- estimateGLMCommonDisp(y,design)
	y <- estimateGLMTrendedDisp(y,design)
	y <- estimateGLMTagwiseDisp(y,design)
	y<-glmFit(y, design)
	yLRT<-glmLRT(y)
	edgeR <- topTags( yLRT,nrow(count) )
#	length(which(edgeR$table$FDR<=0.05))
#	y<-glmLRT(y,contrast=c(-1,1,0,0,0))
#	edgeR <- topTags( y,nrow(count) )
	return(edgeR$table)
}

combined_result2<-function(DESeq,edgeR,bayseq,rawCount,comparGroups=c(0,0,0,1,1,1),rankMethod1=1) {
	DESeq<-as.data.frame(DESeq)
	temp<-row.names(DESeq)
	temp1<-apply(rawCount,1,function(x) log2(mean(x[which(comparGroups==1)])/mean(x[which(comparGroups==0)])))
	temp1<-temp1[temp]
	if (NCOL(DESeq)==6) {
		DESeqResultIndex<-c(2,5,6)
	} else {
		DESeqResultIndex<-c(2,4,5)
	}
	compar1Result<-cbind(DESeq[temp,DESeqResultIndex],edgeR[temp,c(1,4,5)],rawLogFold=temp1,bayseq[temp,])
	rankCompar1Result1<-apply(apply(compar1Result[,c(2,5,8)],2,rank),1,sum)
	if (rankMethod1==1) {
		result<-cbind(compar1Result,apply(compar1Result[,c(2,5,8)],2,rank),rankMethod1=rankCompar1Result1)
	} else {
		rankCompar1Result2<-matrix(NA,nrow=nrow(compar1Result),ncol=(ncol(compar1Result)/2))
		for (x in 1:(ncol(compar1Result)/3)) {
			temp<-sign(compar1Result[,x*3-2])
			temp[temp==0]<-1
			temp<-temp*compar1Result[,x*3-1]
			rankCompar1Result2[,x]<-rank(temp)
		}
		temp<-rankCompar1Result2
		rankCompar1Result2<-apply(rankCompar1Result2,1,sum)
		names(rankCompar1Result2)<-row.names(compar1Result)
		result<-cbind(compar1Result,temp,rankMethod2=rankCompar1Result2)
	}
	colnames(result)[1:12]<-c("log2FoldChange(DESeq2)","pValue(DESeq2)","pAdj(DESeq2)",
	                          "log2FoldChange(edgeR)","pValue(edgeR)","pAdj(edgeR)",
	                          "log2FoldChange(raw)","OneMinusLikelihood(baySeq)","pAdj(baySeq)",
	                          "rank(DESeq)","rank(edgeR)","rank(baySeq)")
	return(result)
}

panel.cor <- function(x, y, digits=2, prefix="", cex.cor) 
{
	usr <- par("usr"); on.exit(par(usr)) 
	par(usr = c(0, 1, 0, 1)) 
	r <- abs(cor(x, y,use="pa",method="sp"))
	txt <- format(c(r, 0.123456789), digits=digits)[1] 
	txt <- paste(prefix, txt, sep="") 
	if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
	
	temp<-options("warn")
	options(warn=-1)
	test <- cor.test(x,y,use="pa",method="sp")
	options(temp)
	
	# borrowed from printCoefmat
	Signif <- symnum(test$p.value, corr = FALSE, na = FALSE, 
			cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
			symbols = c("***", "**", "*", ".", " ")) 
	
	if (cex.cor * r<=1.5) {r<-1.5/cex.cor}
#	r[r<=0.3]<-0.3 #to make the too small r visable
	text(0.5, 0.5, txt, cex = cex.cor * r) 
	text(.8, .8, Signif, cex=cex.cor, col=2) 
}

panel.smooth <-function (x, y, col = par("col"), bg = NA, pch = 16, 
		cex = 1, col.smooth = "red", span = 2/3, iter = 3, ...) 
{
	points(x, y, pch = pch, col = col, bg = bg, cex = cex)
	ok <- is.finite(x) & is.finite(y)
	if (any(ok)) {
		lines(stats::lowess(x[ok], y[ok], f = span, iter = iter), 
				col = col.smooth, ...)
	}
}

venn.diagram1<-function (x, filename, height = 3000, width = 3000, resolution = 500, 
		units = "px", compression = "lzw", na = "stop", main = NULL, 
		sub = NULL, main.pos = c(0.5, 1.05), main.fontface = "plain", 
		main.fontfamily = "serif", main.col = "black", main.cex = 1, 
		main.just = c(0.5, 1), sub.pos = c(0.5, 1.05), sub.fontface = "plain", 
		sub.fontfamily = "serif", sub.col = "black", sub.cex = 1, 
		sub.just = c(0.5, 1), category.names = names(x), force.unique = TRUE, 
		...) 
{
	if (force.unique) {
		for (i in 1:length(x)) {
			x[[i]] <- unique(x[[i]])
		}
	}
	if ("none" == na) {
		x <- x
	}
	else if ("stop" == na) {
		for (i in 1:length(x)) {
			if (any(is.na(x[[i]]))) {
				stop("NAs in dataset", call. = FALSE)
			}
		}
	}
	else if ("remove" == na) {
		for (i in 1:length(x)) {
			x[[i]] <- x[[i]][!is.na(x[[i]])]
		}
	}
	else {
		stop("Invalid na option: valid options are \"none\", \"stop\", and \"remove\"")
	}
	if (0 == length(x) | length(x) > 5) {
		stop("Incorrect number of elements.", call. = FALSE)
	}
	if (1 == length(x)) {
		list.names <- category.names
		if (is.null(list.names)) {
			list.names <- ""
		}
		grob.list <- VennDiagram::draw.single.venn(area = length(x[[1]]), 
				category = list.names, ind = FALSE, ...)
	}
	else if (2 == length(x)) {
		grob.list <- VennDiagram::draw.pairwise.venn(area1 = length(x[[1]]), 
				area2 = length(x[[2]]), cross.area = length(intersect(x[[1]], 
								x[[2]])), category = category.names, ind = FALSE, 
				...)
	}
	else if (3 == length(x)) {
		A <- x[[1]]
		B <- x[[2]]
		C <- x[[3]]
		list.names <- category.names
		nab <- intersect(A, B)
		nbc <- intersect(B, C)
		nac <- intersect(A, C)
		nabc <- intersect(nab, C)
		grob.list <- VennDiagram::draw.triple.venn(area1 = length(A), 
				area2 = length(B), area3 = length(C), n12 = length(nab), 
				n23 = length(nbc), n13 = length(nac), n123 = length(nabc), 
				category = list.names, ind = FALSE, list.order = 1:3, 
				...)
	}
	else if (4 == length(x)) {
		A <- x[[1]]
		B <- x[[2]]
		C <- x[[3]]
		D <- x[[4]]
		list.names <- category.names
		n12 <- intersect(A, B)
		n13 <- intersect(A, C)
		n14 <- intersect(A, D)
		n23 <- intersect(B, C)
		n24 <- intersect(B, D)
		n34 <- intersect(C, D)
		n123 <- intersect(n12, C)
		n124 <- intersect(n12, D)
		n134 <- intersect(n13, D)
		n234 <- intersect(n23, D)
		n1234 <- intersect(n123, D)
		grob.list <- VennDiagram::draw.quad.venn(area1 = length(A), 
				area2 = length(B), area3 = length(C), area4 = length(D), 
				n12 = length(n12), n13 = length(n13), n14 = length(n14), 
				n23 = length(n23), n24 = length(n24), n34 = length(n34), 
				n123 = length(n123), n124 = length(n124), n134 = length(n134), 
				n234 = length(n234), n1234 = length(n1234), category = list.names, 
				ind = FALSE, ...)
	}
	else if (5 == length(x)) {
		A <- x[[1]]
		B <- x[[2]]
		C <- x[[3]]
		D <- x[[4]]
		E <- x[[5]]
		list.names <- category.names
		n12 <- intersect(A, B)
		n13 <- intersect(A, C)
		n14 <- intersect(A, D)
		n15 <- intersect(A, E)
		n23 <- intersect(B, C)
		n24 <- intersect(B, D)
		n25 <- intersect(B, E)
		n34 <- intersect(C, D)
		n35 <- intersect(C, E)
		n45 <- intersect(D, E)
		n123 <- intersect(n12, C)
		n124 <- intersect(n12, D)
		n125 <- intersect(n12, E)
		n134 <- intersect(n13, D)
		n135 <- intersect(n13, E)
		n145 <- intersect(n14, E)
		n234 <- intersect(n23, D)
		n235 <- intersect(n23, E)
		n245 <- intersect(n24, E)
		n345 <- intersect(n34, E)
		n1234 <- intersect(n123, D)
		n1235 <- intersect(n123, E)
		n1245 <- intersect(n124, E)
		n1345 <- intersect(n134, E)
		n2345 <- intersect(n234, E)
		n12345 <- intersect(n1234, E)
		grob.list <- VennDiagram::draw.quintuple.venn(area1 = length(A), 
				area2 = length(B), area3 = length(C), area4 = length(D), 
				area5 = length(E), n12 = length(n12), n13 = length(n13), 
				n14 = length(n14), n15 = length(n15), n23 = length(n23), 
				n24 = length(n24), n25 = length(n25), n34 = length(n34), 
				n35 = length(n35), n45 = length(n45), n123 = length(n123), 
				n124 = length(n124), n125 = length(n125), n134 = length(n134), 
				n135 = length(n135), n145 = length(n145), n234 = length(n234), 
				n235 = length(n235), n245 = length(n245), n345 = length(n345), 
				n1234 = length(n1234), n1235 = length(n1235), n1245 = length(n1245), 
				n1345 = length(n1345), n2345 = length(n2345), n12345 = length(n12345), 
				category = list.names, ind = FALSE, ...)
	}
	else {
		stop("Invalid size of input object")
	}
	if (!is.null(sub)) {
		grob.list <- add.title(gList = grob.list, x = sub, pos = sub.pos, 
				fontface = sub.fontface, fontfamily = sub.fontfamily, 
				col = sub.col, cex = sub.cex)
	}
	if (!is.null(main)) {
		grob.list <- add.title(gList = grob.list, x = main, pos = main.pos, 
				fontface = main.fontface, fontfamily = main.fontfamily, 
				col = main.col, cex = main.cex)
	}
#	if (!is.null(filename)) {
#		current.type <- getOption("bitmapType")
#		if (length(grep("Darwin", Sys.info()["sysname"]))) {
#			options(bitmapType = "quartz")
#		}
#		else {
#			options(bitmapType = "cairo")
#		}
#		tiff(filename = filename, height = height, width = width, 
#				units = units, res = resolution, compression = compression)
	grid.newpage()
	grid.draw(grob.list)
#		dev.off()
#		options(bitmapType = current.type)
	return(1)
#	}
	return(grob.list)
}

plotOverlap<-function(result,cutoffP=0.05,cutoffFC=1,cutoffTopN=100,method=c("p","pf","f","r")) {
#	require(VennDiagram)
	vennResult<-list()
	method<-match.arg(method)
	if (! (method %in% c("pf","p","f","r"))) {
		print("method should be one of pf, p, f, r")
		return(0)
	}
	if (method=="pf") { #p and fold changes cutoff
		
	} else if (method=="p") { #adjust p cutoff
		for (x in 1:3) {
			vennResult[[x]]<-row.names(result)[which(result[,x*3]<=cutoffP)]
			titleText<-paste("Cutoff: adjusted p value <=",cutoffP,sep="")
		}
	} else if (method=="f") { #fold change cutoff
		for (x in 1:3) {
			vennResult[[x]]<-row.names(result)[which(abs(result[,(x*3-2)])>=cutoffFC)]
			titleText<-paste("Cutoff: log2 foldChange >=",cutoffFC,sep="")
		}
	} else if (method=="r") { #top number cutoff
		for (x in 10:12) {
			vennResult[[(x-9)]]<-row.names(result)[which(result[,x]<=cutoffTopN)]
			titleText<-paste("Cutoff: Rank <=",cutoffTopN,sep="")
		}
	}
	names(vennResult)<-c("DESeq2","edgeR","baySeq")
	venn.diagram1(x=vennResult, fill=c("red", "green", "blue"),alpha = c(0.5, 0.5, 0.5),main=titleText,cex=1.5,cat.cex=1.5,main.cex=1.5)
}

heatmap3<-function (x, Rowv = NULL, Colv = if (symm) "Rowv" else NULL, 
		distfun = function(x) as.dist(1 - cor(t(x),use="pa")),balanceColor=F, ColSideLabs,RowSideLabs,showColDendro=T,showRowDendro=T,col=colorRampPalette(c("navy", "white", "firebrick3"))(1024),legendfun,method="complete",ColAxisColors=0,RowAxisColors=0, hclustfun = hclust, reorderfun = function(d, 
				w) reorder(d, w), add.expr,symm = FALSE, revC = identical(Colv, 
				"Rowv"), scale = c("row", "column", "none"), na.rm = TRUE, 
		margins = c(5, 5), ColSideColors, RowSideColors, cexRow = 0.2 + 
				1/log10(nr), cexCol = 0.2 + 1/log10(nc), labRow = NULL, 
		labCol = NULL, main = NULL, xlab = NULL, ylab = NULL, keep.dendro = FALSE, 
		verbose = getOption("verbose"),...) 
{
	scale <- if (missing(scale)) 
				"none"
			else match.arg(scale)
	if (is.data.frame(x)) {x<-as.matrix(x)}
	if (!missing(ColSideColors)) {
		if (is.vector(ColSideColors)) {
			ColSideColors<-cbind(ColSideColors)
		}
	}
	if (!missing(RowSideColors)) {
		if (is.vector(RowSideColors)) {
			RowSideColors<-cbind(RowSideColors)
		}
	}
	if (length(di <- dim(x)) != 2 || !is.numeric(x)) 
		stop("'x' must be a numeric matrix")
	nr <- di[1L]
	nc <- di[2L]
	if (nr <= 1 || nc <= 1) 
		stop("'x' must have at least 2 rows and 2 columns")
	if (!is.numeric(margins) || length(margins) != 2L) 
		stop("'margins' must be a numeric vector of length 2")
	doRdend <- !identical(Rowv, NA)
	doCdend <- !identical(Colv, NA)
	if (!doRdend && identical(Colv, "Rowv")) 
		doCdend <- FALSE
	if (is.null(Rowv)) 
		Rowv <- rowMeans(x, na.rm = na.rm)
	if (is.null(Colv)) 
		Colv <- colMeans(x, na.rm = na.rm)
	if (doRdend) {
		if (inherits(Rowv, "dendrogram")) 
			ddr <- Rowv
		else {
			hcr <- hclustfun(distfun(x),method=method)
			ddr <- as.dendrogram(hcr)
			if (!is.logical(Rowv) || Rowv) 
				ddr <- reorderfun(ddr, Rowv)
		}
		if (nr != length(rowInd <- order.dendrogram(ddr))) 
			stop("row dendrogram ordering gave index of wrong length")
	}
	else rowInd <- 1L:nr
	if (doCdend) {
		if (inherits(Colv, "dendrogram")) 
			ddc <- Colv
		else if (identical(Colv, "Rowv")) {
			if (nr != nc) 
				stop("Colv = \"Rowv\" but nrow(x) != ncol(x)")
			ddc <- ddr
		}
		else {
			hcc <- hclustfun(distfun(if (symm) 
										x
									else t(x)))
			ddc <- as.dendrogram(hcc)
			if (!is.logical(Colv) || Colv) 
				ddc <- reorderfun(ddc, Colv)
		}
		if (nc != length(colInd <- order.dendrogram(ddc))) 
			stop("column dendrogram ordering gave index of wrong length")
	}
	else colInd <- 1L:nc
	x <- x[rowInd, colInd]
	labRow <- if (is.null(labRow)) 
				if (is.null(rownames(x))) 
					(1L:nr)[rowInd]
				else rownames(x)
			else labRow[rowInd]
	labCol <- if (is.null(labCol)) 
				if (is.null(colnames(x))) 
					(1L:nc)[colInd]
				else colnames(x)
			else labCol[colInd]
	if (scale == "row") {
		x <- sweep(x, 1L, rowMeans(x, na.rm = na.rm), check.margin = FALSE)
		sx <- apply(x, 1L, sd, na.rm = na.rm)
		x <- sweep(x, 1L, sx, "/", check.margin = FALSE)
	}
	else if (scale == "column") {
		x <- sweep(x, 2L, colMeans(x, na.rm = na.rm), check.margin = FALSE)
		sx <- apply(x, 2L, sd, na.rm = na.rm)
		x <- sweep(x, 2L, sx, "/", check.margin = FALSE)
	}
	lmat <- rbind(c(NA, 3), 2:1)
#	lwid <- c(if (doRdend) 1 else 0.05, 4)
	lwid <- c(1, 4)
#	lhei <- c((if (doCdend) 1 else 0.05) + if (!is.null(main)) 0.2 else 0, 
#			4)
	lhei <- c( 1 + if (!is.null(main)) 0.2 else 0,4)
	if (!missing(ColSideColors)) {
		if (!is.character(ColSideColors) || nrow(ColSideColors) != 
				nc) 
			stop("'ColSideColors' must be a character vector or matrix of length ncol(x)")
		lmat <- rbind(lmat[1, ] + 1, c(NA, 1), lmat[2, ] + 1)
		lhei <- c(lhei[1L], 0.2, lhei[2L])
	}
	if (!missing(RowSideColors)) {
		if (!is.character(RowSideColors) || nrow(RowSideColors) != 
				nr) 
			stop("'RowSideColors' must be a character vector or matrix of length nrow(x)")
		lmat <- cbind(lmat[, 1] + 1, c(rep(NA, nrow(lmat) - 1), 
						1), lmat[, 2] + 1)
		lwid <- c(lwid[1L], 0.2, lwid[2L])
	}
	lmat<-lmat+1
	lmat[is.na(lmat)] <- 0
	lmat[1,1]<-1
	if (verbose) {
		cat("layout: widths = ", lwid, ", heights = ", lhei, 
				"; lmat=\n")
		print(lmat)
	}
	dev.hold()
	on.exit(dev.flush())
	op <- par(no.readonly = TRUE)
	on.exit(par(op), add = TRUE)
	
	#balanceColor
	if (balanceColor) {
		if (abs(max(x,na.rm=T))>=abs(min(x,na.rm=T))) {
			cut.off<-round(quantile(1:length(col),probs=1-(abs(max(x,na.rm=T))+abs(min(x,na.rm=T)))/(2*abs(max(x,na.rm=T)))))
			col<-col[cut.off:length(col)]
		} else {
			cut.off<-round(quantile(1:length(col),probs=(abs(max(x,na.rm=T))+abs(min(x,na.rm=T)))/(2*abs(min(x,na.rm=T)))))
			col<-col[1:cut.off]
		}
	}
	
	layout(lmat, widths = lwid, heights = lhei, respect = TRUE)
	if (!missing(legendfun)) {
		par(mar = c(0, 0, 0, 0))
		legendfun()
	} else {
		par(mar = c(5, 1, 1, 0))
		dummy.x <- seq(min(x, na.rm = TRUE), max(x, na.rm = TRUE), 
				length = length(col))
		dummy.z <- matrix(dummy.x, ncol = 1)
		image(x = dummy.x, y = 1, z = dummy.z, yaxt = "n", 
				col = col,cex.axis=cexCol,xlab="")
	}
	if (!missing(RowSideColors)) {
		par(mar = c(margins[1L], 0, 0, 0.5))
		if (revC) {
			rsc = RowSideColors[rev(rowInd),,drop=F]
		} else {
			rsc = RowSideColors[rowInd,,drop=F]
		}
		rsc.colors = matrix()
		rsc.names = names(table(rsc))
		rsc.i = 1
		for (rsc.name in rsc.names) {
			rsc.colors[rsc.i] = rsc.name
			rsc[rsc == rsc.name] = rsc.i
			rsc.i = rsc.i + 1
		}
		rsc = matrix(as.numeric(rsc), nrow = dim(rsc)[1])
		image(t(rsc), col = as.vector(rsc.colors), axes = FALSE)
		
		if (missing(RowSideLabs)) {
			RowSideLabs<-colnames(RowSideColors)
		}
		if (dim(rsc)[2]==1) {
			axis(1, 0, RowSideLabs, las = 2, tick = FALSE)
		} else {
			axis(1, 0:(dim(rsc)[2] - 1)/(dim(rsc)[2] - 1), RowSideLabs, 
					las = 2, tick = FALSE)
		}
	}
	if (!missing(ColSideColors)) {
		par(mar = c(0.5, 0, 0, margins[2L]))
		csc = ColSideColors[colInd,,drop=F]
		csc.colors = matrix()
		csc.names = names(table(csc))
		csc.i = 1
		for (csc.name in csc.names) {
			csc.colors[csc.i] = csc.name
			csc[csc == csc.name] = csc.i
			csc.i = csc.i + 1
		}
		csc = matrix(as.numeric(csc), nrow = dim(csc)[1])
		image(csc, col = as.vector(csc.colors), axes = FALSE)
		if (missing(ColSideLabs)) {
			ColSideLabs<-colnames(ColSideColors)
		}
		if (dim(csc)[2]==1) {
			axis(4, 0, ColSideLabs,las = 2, tick = FALSE)
		} else {
			axis(4, 0:(dim(csc)[2] - 1)/(dim(csc)[2] - 1), ColSideLabs,las = 2, tick = FALSE)
		}
	}
	par(mar = c(margins[1L], 0, 0, margins[2L]))
	if (!symm || scale != "none") 
		x <- t(x)
	if (revC) {
		iy <- nr:1
		if (doRdend)
			ddr <- rev(ddr)
		x <- x[, iy]
	}
	else iy <- 1L:nr
	image(1L:nc, 1L:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 + 
					c(0, nr), axes = FALSE, xlab = "", ylab = "", col=col,...)
	if (!missing(ColSideColors) & ColAxisColors!=0) {
		mtext(1, at=1L:nc, text = labCol, las = 2, line = 0.5,cex = cexCol,col=ColSideColors[colInd,ColAxisColors])
	} else {
		axis(1, 1L:nc, labels = labCol, las = 2, line = -0.5, tick = 0,cex.axis = cexCol)
	}
	
	if (!is.null(xlab)) 
		mtext(xlab, side = 1, line = margins[1L] - 1.25)
	if (!missing(RowSideColors) & RowAxisColors!=0) {
		mtext(4, at=iy, text = labRow, las = 2, line = 0.5,cex = cexRow,col=RowSideColors[rowInd,RowAxisColors])
	} else {
		axis(4, iy, labels = labRow, las = 2, line = -0.5, tick = 0, cex.axis = cexRow)
	}
	if (!is.null(ylab)) 
		mtext(ylab, side = 4, line = margins[2L] - 1.25)
	if (!missing(add.expr)) 
		eval(substitute(add.expr))
	par(mar = c(margins[1L], 0, 0, 0))
	if (doRdend & showRowDendro) 
		plot(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
	else frame()
	par(mar = c(0, 0, if (!is.null(main)) 1 else 0, margins[2L]))
	if (doCdend & showColDendro) {
		plot(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
	}	else if (!is.null(main)) 
		frame()
	if (!is.null(main)) {
		par(xpd = NA)
		title(main, cex.main = 1.5 * op[["cex.main"]])
	}
	invisible(list(rowInd = rowInd, colInd = colInd, Rowv = if (keep.dendro && 
							doRdend) ddr, Colv = if (keep.dendro && doCdend) ddc))
}

history1<-function (max.show = 25, reverse = FALSE, pattern, ...) 
{
	file1 <- tempfile("Rrawhist")
	savehistory(file1)
	rawhist <- readLines(file1)
	unlink(file1)
	if (!missing(pattern)) 
		rawhist <- unique(grep(pattern, rawhist, value = TRUE, 
						...))
	inds<-length(rawhist)
	return(rawhist[inds])
}

density.plot<-function(data,width=c(-1,15),col="",do.legend=T,log.density=F,main="",adjust=1,cex.main=1,cex.axis=1,legend.text="",cex.lab=1,cex.legend=1,legend.x=0.5,lwd=2,ylab="Density",xlab="",legend.lwd=2,...) {
	if (log.density==T) {data<-log2(data)}
	if (legend.text[1]=="") {legend.text=colnames(data)}
	sample<-length(data[1,])
	temp<-0
	for (x in 1:sample) {temp[x]<-max(density(data[,x],na.rm=T)$y)}
	temp<-max(temp)
	if (col[1]=="") {col=rainbow(sample)}
	plot(width,c(0,temp),type="n",ylab=ylab,xlab=xlab,main=main,cex.axis=cex.axis,cex.main=cex.main,cex.lab=cex.lab,las=1,...)
	for (x in 1:sample) {lines(density(data[,x],adjust=adjust,na.rm=T),lwd=lwd,col=rep(col,round(sample/length(col))+1)[x])}
	if (do.legend==T) {
		legend(width[2]*legend.x,temp,legend=legend.text,lwd=legend.lwd,col=col,bty = "n",cex=cex.legend)
	}
}