--- a
+++ b/modas/utils/CMplot.r
@@ -0,0 +1,1631 @@
+#Version:3.3.3
+#Data: 2018/12/11
+#Author: Lilin Yin
+
+CMplot <- function(
+	Pmap,
+	col=c("#377EB8", "#4DAF4A", "#984EA3", "#FF7F00"),
+	bin.size=1e6,
+	bin.max=NULL,
+	pch=19,
+	band=1,
+	cir.band=0.5,
+	H=1.5,
+	ylim=NULL,
+	cex.axis=1,
+	plot.type="b",
+	multracks=FALSE,
+	cex=c(0.5,1,1),
+	r=0.3,
+	xlab="Chromosome",
+	ylab=expression(-log[10](italic(p))),
+	xaxs="i",
+	yaxs="r",
+	outward=FALSE,
+	threshold = NULL, 
+	threshold.col="red",
+	threshold.lwd=1,
+	threshold.lty=2,
+	amplify= TRUE,
+	chr.labels=NULL,
+	signal.cex = 1.5,
+	signal.pch = 19,
+	signal.col="red",
+	signal.line=1,
+	cir.chr=TRUE,
+	cir.chr.h=1.5,
+	chr.den.col=c("darkgreen", "yellow", "red"),
+	cir.legend=TRUE,
+	cir.legend.cex=0.6,
+	cir.legend.col="black",
+	LOG10=TRUE,
+	box=FALSE,
+	conf.int.col="grey",
+	file.output=TRUE,
+	file="jpg",
+	dpi=300,
+	memo="",
+	verbose=TRUE
+)
+{	
+
+	#plot a circle with a radius of r
+	circle.plot <- function(myr,type="l",x=NULL,lty=1,lwd=1,col="black",add=TRUE,n.point=1000)
+	{
+		curve(sqrt(myr^2-x^2),xlim=c(-myr,myr),n=n.point,ylim=c(-myr,myr),type=type,lty=lty,col=col,lwd=lwd,add=add)
+		curve(-sqrt(myr^2-x^2),xlim=c(-myr,myr),n=n.point,ylim=c(-myr,myr),type=type,lty=lty,col=col,lwd=lwd,add=TRUE)
+	}
+	
+	Densitplot <- function(
+		map,
+		col=c("darkgreen", "yellow", "red"),
+		main="SNP Density",
+		bin=1e6,
+		band=3,
+		width=5,
+		legend.len=10,
+		legend.max=NULL,
+		legend.pt.cex=3,
+		legend.cex=1,
+		legend.y.intersp=1,
+		legend.x.intersp=1,
+		plot=TRUE
+	)
+	{
+		map <- as.matrix(map)
+		map <- map[!is.na(map[, 2]), ]
+		map <- map[!is.na(as.numeric(map[, 3])), ]
+		map <- map[map[, 2] != 0, ]
+		#map <- map[map[, 3] != 0, ]
+		options(warn = -1)
+		max.chr <- max(as.numeric(map[, 2]), na.rm=TRUE)
+		if(is.infinite(max.chr))	max.chr <- 0
+		map.xy.index <- which(!as.numeric(map[, 2]) %in% c(0 : max.chr))
+		if(length(map.xy.index) != 0){
+			chr.xy <- unique(map[map.xy.index, 2])
+			for(i in 1:length(chr.xy)){
+				map[map[, 2] == chr.xy[i], 2] <- max.chr + i
+			}
+		}
+		map <- map[order(as.numeric(map[, 2]), as.numeric(map[, 3])), ]
+		chr <- as.numeric(map[, 2])
+		pos <- as.numeric(map[, 3])
+		chr.num <- unique(chr)
+		chorm.maxlen <- max(pos)
+		if(plot)	plot(NULL, xlim=c(0, chorm.maxlen + chorm.maxlen/10), ylim=c(0, length(chr.num) * band + band), main=main,axes=FALSE, xlab="", ylab="", xaxs="i", yaxs="i")
+		pos.x <- list()
+		col.index <- list()
+		maxbin.num <- NULL
+		for(i in 1 : length(chr.num)){
+			pos.x[[i]] <- pos[which(chr == chr.num[i])]
+			cut.len <- ceiling((max(pos.x[[i]]) - min(pos.x[[i]])) / bin)
+			if(cut.len <= 1){
+				maxbin.num <- c(maxbin.num,length(pos.x[[i]]))
+               	col.index[[i]] <- rep(length(pos.x[[i]]), length(pos.x[[i]]))
+			}else{
+				cut.r <- cut(pos.x[[i]], cut.len, labels=FALSE)
+				eachbin.num <- table(cut.r)
+				maxbin.num <- c(maxbin.num, max(eachbin.num))
+				col.index[[i]] <- rep(eachbin.num, eachbin.num)
+			}
+		}
+		Maxbin.num <- max(maxbin.num)
+		maxbin.num <- Maxbin.num
+		if(!is.null(legend.max)){
+			maxbin.num <- legend.max
+		}
+		col=colorRampPalette(col)(maxbin.num)
+		col.seg=NULL
+		for(i in 1 : length(chr.num)){
+			if(plot)	polygon(c(0, 0, max(pos.x[[i]]), max(pos.x[[i]])), 
+				c(-width/5 - band * (i - length(chr.num) - 1), width/5 - band * (i - length(chr.num) - 1), 
+				width/5 - band * (i - length(chr.num) - 1), -width/5 - band * (i - length(chr.num) - 1)), col="grey", border="grey")
+			if(!is.null(legend.max)){
+				if(legend.max < Maxbin.num){
+					col.index[[i]][col.index[[i]] > legend.max] <- legend.max
+				}
+			}
+			col.seg <- c(col.seg, col[round(col.index[[i]] * length(col) / maxbin.num)])
+			if(plot)	segments(pos.x[[i]], -width/5 - band * (i - length(chr.num) - 1), pos.x[[i]], width/5 - band * (i - length(chr.num) - 1), 
+			col=col[round(col.index[[i]] * length(col) / maxbin.num)], lwd=1)
+		}
+		if(length(map.xy.index) != 0){
+			for(i in 1:length(chr.xy)){
+				chr.num[chr.num == max.chr + i] <- chr.xy[i]
+			}
+		}
+		chr.num <- rev(chr.num)
+		if(plot){
+			if(max.chr == 0)	mtext(at=seq(band, length(chr.num) * band, band),text=chr.num, side=2, las=2, font=1, cex=cex.axis*0.6, line=0.2)
+			if(max.chr != 0)	mtext(at=seq(band, length(chr.num) * band, band),text=paste("Chr", chr.num, sep=""), side=2, las=2, font=1, cex=cex.axis*0.6, line=0.2)
+		}
+		if(plot)	axis(3, at=seq(0, chorm.maxlen, length=10), labels=paste(round((seq(0, chorm.maxlen, length=10)) / 1e6, 0), "Mb", sep=""),
+			font=1, cex.axis=cex.axis*0.8, tck=0.01, lwd=2, padj=1.2)
+		# image(c(chorm.maxlen-chorm.maxlen * legend.width / 20 , chorm.maxlen), 
+		# round(seq(band - width/5, (length(chr.num) * band + band) * legend.height / 2 , length=maxbin.num+1), 2), 
+		# t(matrix(0 : maxbin.num)), col=c("white", rev(heat.colors(maxbin.num))), add=TRUE)
+				
+		if(maxbin.num <= legend.len)	legend.len <- maxbin.num		
+		
+		legend.y <- round(seq(0, maxbin.num, length=legend.len))
+		len <- legend.y[2]
+		legend.y <- seq(0, maxbin.num, len)
+		if(!is.null(legend.max)){
+			if(legend.max < Maxbin.num){
+				if(!maxbin.num %in% legend.y){
+					legend.y <- c(legend.y, paste(">=", maxbin.num, sep=""))
+					legend.y.col <- c(legend.y[c(-1, -length(legend.y))], maxbin.num)
+				}else{
+					legend.y[length(legend.y)] <- paste(">=", maxbin.num, sep="")
+					legend.y.col <- c(legend.y[c(-1, -length(legend.y))], maxbin.num)
+				}
+			}else{
+				if(!maxbin.num %in% legend.y){
+					legend.y <- c(legend.y, maxbin.num)
+				}
+				legend.y.col <- c(legend.y[-1])
+			}
+		}else{
+			if(!maxbin.num %in% legend.y){
+				legend.y <- c(legend.y, paste(">", max(legend.y), sep=""))
+				legend.y.col <- c(legend.y[c(-1, -length(legend.y))], maxbin.num)
+			}else{
+				legend.y.col <- c(legend.y[-1])
+			}
+		}
+		legend.y.col <- as.numeric(legend.y.col)
+		legend.col <- c("grey", col[round(legend.y.col * length(col) / maxbin.num)])
+		if(plot)	legend(x=(chorm.maxlen + chorm.maxlen/100), y=( -width/2.5 - band * (length(chr.num) - length(chr.num) - 1)), title="", legend=legend.y, pch=15, pt.cex = legend.pt.cex, col=legend.col,
+			cex=legend.cex, bty="n", y.intersp=legend.y.intersp, x.intersp=legend.x.intersp, yjust=0, xjust=0, xpd=TRUE)
+		if(!plot)	return(list(den.col=col.seg, legend.col=legend.col, legend.y=legend.y))
+	}
+
+	if(sum(plot.type %in% "b")==1) plot.type=c("c","m","q","d")
+
+	taxa=colnames(Pmap)[-c(1:3)]
+	if(!is.null(memo) && memo != "")	memo <- paste("_", memo, sep="")
+	if(length(taxa) == 0)	taxa <- "Index"
+	taxa <- paste(taxa, memo, sep="")
+
+	#SNP-Density plot
+	if("d" %in% plot.type){
+		if(verbose)	print("SNP_Density Plotting...")
+		if(file.output){
+			if(file=="jpg")	jpeg(paste("SNP_Density.",paste(taxa,collapse="."),".jpg",sep=""), width = 9*dpi,height=7*dpi,res=dpi,quality = 100)
+			if(file=="pdf")	pdf(paste("SNP_Density.",paste(taxa,collapse="."),".pdf",sep=""), width = 9,height=7)
+			if(file=="tiff")	tiff(paste("SNP_Density.",paste(taxa,collapse="."),".tiff",sep=""), width = 9*dpi,height=7*dpi,res=dpi)
+			par(xpd=TRUE)
+		}else{
+			if(is.null(dev.list()))	dev.new(width = 9,height=7)
+			par(xpd=TRUE)
+		}
+
+		Densitplot(map=Pmap[,c(1:3)], col=col, bin=bin.size, legend.max=bin.max, main=paste("The number of SNPs within ", bin.size/1e6, "Mb window size", sep=""))
+		if(file.output)	dev.off()
+	}
+
+	if(length(plot.type) !=1 | (!"d" %in% plot.type)){
+	
+		#order Pmap by the name of SNP
+		#Pmap=Pmap[order(Pmap[,1]),]
+		Pmap <- as.matrix(Pmap)
+
+		#delete the column of SNPs names
+		Pmap <- Pmap[,-1]
+		
+		#scale and adjust the parameters
+		cir.chr.h <- cir.chr.h/5
+		cir.band <- cir.band/5
+		if(!is.null(threshold)){
+			threshold.col <- rep(threshold.col,length(threshold))
+			threshold.lwd <- rep(threshold.lwd,length(threshold))
+			threshold.lty <- rep(threshold.lty,length(threshold))
+			signal.col <- rep(signal.col,length(threshold))
+			signal.pch <- rep(signal.pch,length(threshold))
+			signal.cex <- rep(signal.cex,length(threshold))
+		}
+		if(length(cex)!=3) cex <- rep(cex,3)
+		if(!is.null(ylim)){
+			if(length(ylim)==1) ylim <- c(0,ylim)
+		}
+		
+		if(is.null(conf.int.col))	conf.int.col <- NA
+		if(is.na(conf.int.col)){
+			conf.int=FALSE
+		}else{
+			conf.int=TRUE
+		}
+
+		#get the number of traits
+		R=ncol(Pmap)-2
+
+		#replace the non-euchromosome
+		options(warn = -1)
+		numeric.chr <- as.numeric(Pmap[, 1])
+		options(warn = 0)
+		max.chr <- max(numeric.chr, na.rm=TRUE)
+		if(is.infinite(max.chr))	max.chr <- 0
+		map.xy.index <- which(!numeric.chr %in% c(0:max.chr))
+		if(length(map.xy.index) != 0){
+			chr.xy <- unique(Pmap[map.xy.index, 1])
+			for(i in 1:length(chr.xy)){
+				Pmap[Pmap[, 1] == chr.xy[i], 1] <- max.chr + i
+			}
+		}
+
+		Pmap <- matrix(as.numeric(Pmap), nrow(Pmap))
+		Pmap[is.na(Pmap)] <- 0
+
+		#order the GWAS results by chromosome and position
+		Pmap <- Pmap[order(Pmap[, 1], Pmap[,2]), ]
+
+		#get the index of chromosome
+		chr <- unique(Pmap[,1])
+		chr.ori <- chr
+		if(length(map.xy.index) != 0){
+			for(i in 1:length(chr.xy)){
+				chr.ori[chr.ori == max.chr + i] <- chr.xy[i]
+			}
+		}
+
+		pvalueT <- as.matrix(Pmap[,-c(1:2)])
+		pvalue.pos <- Pmap[, 2]
+		p0.index <- Pmap[, 1] == 0
+		if(sum(p0.index) != 0){
+			pvalue.pos[p0.index] <- 1:sum(p0.index)
+		}
+		pvalue.pos.list <- tapply(pvalue.pos, Pmap[, 1], list)
+		
+		#scale the space parameter between chromosomes
+		if(!missing(band)){
+			band <- floor(band*(sum(sapply(pvalue.pos.list, max))/100))
+		}else{
+			band <- floor((sum(sapply(pvalue.pos.list, max))/100))
+		}
+		if(band==0)	band=1
+		
+		if(LOG10){
+			pvalueT[pvalueT <= 0] <- 1
+			pvalueT[pvalueT > 1] <- 1
+		}
+
+		#set the colors for the plot
+		#palette(heat.colors(1024)) #(heatmap)
+		#T=floor(1024/max(pvalue))
+		#plot(pvalue,pch=19,cex=0.6,col=(1024-floor(pvalue*T)))
+		if(is.vector(col)){
+			col <- matrix(col,R,length(col),byrow=TRUE)
+		}
+		if(is.matrix(col)){
+			#try to transform the colors into matrix for all traits
+			col <- matrix(as.vector(t(col)),R,dim(col)[2],byrow=TRUE)
+		}
+
+		Num <- as.numeric(table(Pmap[,1]))
+		Nchr <- length(Num)
+		N <- NULL
+
+		#set the colors for each traits
+		for(i in 1:R){
+			colx <- col[i,]
+			colx <- colx[!is.na(colx)]
+			N[i] <- ceiling(Nchr/length(colx))
+		}
+		
+		#insert the space into chromosomes and return the midpoint of each chromosome
+		ticks <- NULL
+		pvalue.posN <- NULL
+		#pvalue <- pvalueT[,j]
+		for(i in 0:(Nchr-1)){
+			if (i==0){
+				#pvalue <- append(pvalue,rep(Inf,band),after=0)
+				pvalue.posN <- pvalue.pos.list[[i+1]] + band
+				ticks[i+1] <- max(pvalue.posN)-floor(max(pvalue.pos.list[[i+1]])/2)
+			}else{
+				#pvalue <- append(pvalue,rep(Inf,band),after=sum(Num[1:i])+i*band)
+				pvalue.posN <- c(pvalue.posN, max(pvalue.posN) + band + pvalue.pos.list[[i+1]])
+				ticks[i+1] <- max(pvalue.posN)-floor(max(pvalue.pos.list[[i+1]])/2)
+			}
+		}
+		pvalue.posN.list <- tapply(pvalue.posN, Pmap[, 1], list)
+		#NewP[[j]] <- pvalue
+		
+		#merge the pvalues of traits by column
+		if(LOG10){
+			logpvalueT <- -log10(pvalueT)
+		}else{
+			logpvalueT <- pvalueT
+		}
+
+		add <- list()
+		for(i in 1:R){
+			colx <- col[i,]
+			colx <- colx[!is.na(colx)]
+			add[[i]] <- c(Num,rep(0,N[i]*length(colx)-Nchr))
+		}
+
+		TotalN <- max(pvalue.posN)
+
+		if(length(chr.den.col) > 1){
+			cir.density=TRUE
+			den.fold <- 20
+			density.list <- Densitplot(map=Pmap[,c(1,1,2)], col=chr.den.col, plot=FALSE, bin=bin.size, legend.max=bin.max)
+			#list(den.col=col.seg, legend.col=legend.col, legend.y=legend.y)
+		}else{
+			cir.density=FALSE
+		}
+		
+		signal.line.index <- NULL
+		if(!is.null(threshold)){
+			if(!is.null(signal.line)){
+				for(l in 1:R){
+					if(LOG10){
+						signal.line.index <- c(signal.line.index,which(pvalueT[,l] < min(threshold)))
+					}else{
+						signal.line.index <- c(signal.line.index,which(pvalueT[,l] > max(threshold)))
+					}
+				}
+				signal.line.index <- unique(signal.line.index)
+			}
+		}
+		signal.line.index <- pvalue.posN[signal.line.index]
+	}
+	
+	#plot circle Manhattan
+	if("c" %in% plot.type){
+		#print("Starting Circular-Manhattan plot!",quote=F)
+		if(file.output){
+			if(file=="jpg")	jpeg(paste("Circular-Manhattan.",paste(taxa,collapse="."),".jpg",sep=""), width = 8*dpi,height=8*dpi,res=dpi,quality = 100)
+			if(file=="pdf")	pdf(paste("Circular-Manhattan.",paste(taxa,collapse="."),".pdf",sep=""), width = 10,height=10)
+			if(file=="tiff")	tiff(paste("Circular-Manhattan.",paste(taxa,collapse="."),".tiff",sep=""), width = 8*dpi,height=8*dpi,res=dpi)
+			par(pty="s", xpd=TRUE, mar=c(1,1,1,1))
+		}
+		if(!file.output){
+			if(is.null(dev.list()))	dev.new(width=8, height=8)
+			par(pty="s", xpd=TRUE)
+		}
+		RR <- r+H*R+cir.band*R
+		if(cir.density){
+			plot(NULL,xlim=c(1.05*(-RR-4*cir.chr.h),1.1*(RR+4*cir.chr.h)),ylim=c(1.05*(-RR-4*cir.chr.h),1.1*(RR+4*cir.chr.h)),axes=FALSE,xlab="",ylab="")
+		}else{
+			plot(NULL,xlim=c(1.05*(-RR-4*cir.chr.h),1.05*(RR+4*cir.chr.h)),ylim=c(1.05*(-RR-4*cir.chr.h),1.05*(RR+4*cir.chr.h)),axes=FALSE,xlab="",ylab="")
+		}
+		if(!is.null(signal.line)){
+			if(!is.null(signal.line.index)){
+				X1chr <- (RR)*sin(2*pi*(signal.line.index-round(band/2))/TotalN)
+				Y1chr <- (RR)*cos(2*pi*(signal.line.index-round(band/2))/TotalN)
+				X2chr <- (r)*sin(2*pi*(signal.line.index-round(band/2))/TotalN)
+				Y2chr <- (r)*cos(2*pi*(signal.line.index-round(band/2))/TotalN)
+				segments(X1chr,Y1chr,X2chr,Y2chr,lty=2,lwd=signal.line,col="grey")
+			}
+		}
+		for(i in 1:R){
+		
+			#get the colors for each trait
+			colx <- col[i,]
+			colx <- colx[!is.na(colx)]
+			
+			#debug
+			#print(colx)
+			
+			if(verbose)	print(paste("Circular_Manhattan Plotting ",taxa[i],"...",sep=""))
+			pvalue <- pvalueT[,i]
+			logpvalue <- logpvalueT[,i]
+			if(is.null(ylim)){
+				if(LOG10){
+					Max <- ceiling(-log10(min(pvalue[pvalue!=0])))
+					Min <- 0
+				}else{
+					Max <- ceiling(max(pvalue[pvalue!=Inf]))
+					if(abs(Max)<=1)	Max <- max(pvalue[pvalue!=Inf])
+					Min <- floor(min(pvalue[pvalue!=Inf]))
+					if(abs(Min)<=1)	Min <- min(pvalue[pvalue!=Inf])
+				}
+			}else{
+				Max <- ylim[2]
+				Min <- ylim[1]
+			}
+			Cpvalue <- (H*(logpvalue-Min))/(Max-Min)
+			if(outward==TRUE){
+				if(cir.chr==TRUE){
+					
+					#plot the boundary which represents the chromosomes
+					polygon.num <- 1000
+					for(k in 1:length(chr)){
+						if(k==1){
+							polygon.index <- seq(round(band/2)+1,-round(band/2)+max(pvalue.posN.list[[1]]), length=polygon.num)
+							#change the axis from right angle into circle format
+							X1chr=(RR)*sin(2*pi*(polygon.index)/TotalN)
+							Y1chr=(RR)*cos(2*pi*(polygon.index)/TotalN)
+							X2chr=(RR+cir.chr.h)*sin(2*pi*(polygon.index)/TotalN)
+							Y2chr=(RR+cir.chr.h)*cos(2*pi*(polygon.index)/TotalN)
+							if(is.null(chr.den.col)){
+								polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=rep(colx,ceiling(length(chr)/length(colx)))[k],border=rep(colx,ceiling(length(chr)/length(colx)))[k])	
+							}else{
+								if(cir.density){
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
+								}else{
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
+								}
+							}
+						}else{
+							polygon.index <- seq(1+round(band/2)+max(pvalue.posN.list[[k-1]]),-round(band/2)+max(pvalue.posN.list[[k]]), length=polygon.num)
+							X1chr=(RR)*sin(2*pi*(polygon.index)/TotalN)
+							Y1chr=(RR)*cos(2*pi*(polygon.index)/TotalN)
+							X2chr=(RR+cir.chr.h)*sin(2*pi*(polygon.index)/TotalN)
+							Y2chr=(RR+cir.chr.h)*cos(2*pi*(polygon.index)/TotalN)
+							if(is.null(chr.den.col)){
+								polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=rep(colx,ceiling(length(chr)/length(colx)))[k],border=rep(colx,ceiling(length(chr)/length(colx)))[k])
+							}else{
+								if(cir.density){
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
+								}else{
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
+								}
+							}		
+						}
+					}
+					
+					if(cir.density){
+
+						segments(
+							(RR)*sin(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							(RR)*cos(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							(RR+cir.chr.h)*sin(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							(RR+cir.chr.h)*cos(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							col=density.list$den.col, lwd=0.1
+						)
+						legend(
+							x=RR+4*cir.chr.h,
+							y=(RR+4*cir.chr.h)/2,
+							title="", legend=density.list$legend.y, pch=15, pt.cex = 3, col=density.list$legend.col,
+							cex=1, bty="n",
+							y.intersp=1,
+							x.intersp=1,
+							yjust=0.5, xjust=0, xpd=TRUE
+						)
+						
+					}
+					
+					# XLine=(RR+cir.chr.h)*sin(2*pi*(1:TotalN)/TotalN)
+					# YLine=(RR+cir.chr.h)*cos(2*pi*(1:TotalN)/TotalN)
+					# lines(XLine,YLine,lwd=1.5)
+					if(cir.density){
+						circle.plot(myr=RR+cir.chr.h,lwd=1.5,add=TRUE,col='grey')
+						circle.plot(myr=RR,lwd=1.5,add=TRUE,col='grey')
+					}else{
+						circle.plot(myr=RR+cir.chr.h,lwd=1.5,add=TRUE)
+						circle.plot(myr=RR,lwd=1.5,add=TRUE)
+					}
+
+				}
+				
+				X=(Cpvalue+r+H*(i-1)+cir.band*(i-1))*sin(2*pi*(pvalue.posN-round(band/2))/TotalN)
+				Y=(Cpvalue+r+H*(i-1)+cir.band*(i-1))*cos(2*pi*(pvalue.posN-round(band/2))/TotalN)
+				points(X,Y,pch=19,cex=cex[1],col=rep(rep(colx,N[i]),add[[i]]))
+				
+				#plot the legend for each trait
+				if(cir.legend==TRUE){
+					#try to get the number after radix point
+					if((Max-Min) > 1) {
+						round.n=2
+					}else{
+						round.n=nchar(as.character(10^(-ceiling(-log10(Max)))))-1
+					}
+					segments(0,r+H*(i-1)+cir.band*(i-1),0,r+H*i+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					segments(0,r+H*(i-1)+cir.band*(i-1),H/20,r+H*(i-1)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-1)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0.75)+cir.band*(i-1),H/20,r+H*(i-0.75)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0.75)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0.5)+cir.band*(i-1),H/20,r+H*(i-0.5)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0.5)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0.25)+cir.band*(i-1),H/20,r+H*(i-0.25)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0.25)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0)+cir.band*(i-1),H/20,r+H*(i-0)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					text(-r/15,r+H*(i-0.94)+cir.band*(i-1),round(Min+(Max-Min)*0,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.75)+cir.band*(i-1),round(Min+(Max-Min)*0.25,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.5)+cir.band*(i-1),round(Min+(Max-Min)*0.5,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.25)+cir.band*(i-1),round(Min+(Max-Min)*0.75,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.06)+cir.band*(i-1),round(Min+(Max-Min)*1,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+				}
+				
+				if(!is.null(threshold)){
+					if(sum(threshold!=0)==length(threshold)){
+						for(thr in 1:length(threshold)){
+							significantline1=ifelse(LOG10, H*(-log10(threshold[thr])-Min)/(Max-Min), H*(threshold[thr]-Min)/(Max-Min))
+							#s1X=(significantline1+r+H*(i-1)+cir.band*(i-1))*sin(2*pi*(0:TotalN)/TotalN)
+							#s1Y=(significantline1+r+H*(i-1)+cir.band*(i-1))*cos(2*pi*(0:TotalN)/TotalN)
+							if(significantline1<H){
+								#lines(s1X,s1Y,type="l",col=threshold.col,lwd=threshold.col,lty=threshold.lty)
+								circle.plot(myr=(significantline1+r+H*(i-1)+cir.band*(i-1)),col=threshold.col[thr],lwd=threshold.lwd[thr],lty=threshold.lty[thr])
+							}else{
+								warning(paste("No significant points for ",taxa[i]," pass the threshold level using threshold=",threshold[thr],"!",sep=""))
+							}
+						}
+					}
+				}
+				
+				if(!is.null(threshold)){
+					if(sum(threshold!=0)==length(threshold)){
+						if(amplify==TRUE){
+							if(LOG10){
+								threshold <- sort(threshold)
+								significantline1=H*(-log10(max(threshold))-Min)/(Max-Min)
+							}else{
+								threshold <- sort(threshold, decreasing=TRUE)
+								significantline1=H*(min(threshold)-Min)/(Max-Min)
+							}
+							
+							p_amp.index <- which(Cpvalue>=significantline1)
+							HX1=(Cpvalue[p_amp.index]+r+H*(i-1)+cir.band*(i-1))*sin(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+							HY1=(Cpvalue[p_amp.index]+r+H*(i-1)+cir.band*(i-1))*cos(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+							
+							#cover the points that exceed the threshold with the color "white"
+							points(HX1,HY1,pch=19,cex=cex[1],col="white")
+							
+								for(ll in 1:length(threshold)){
+									if(ll == 1){
+										if(LOG10){
+											significantline1=H*(-log10(threshold[ll])-Min)/(Max-Min)
+										}else{
+											significantline1=H*(threshold[ll]-Min)/(Max-Min)
+										}
+										p_amp.index <- which(Cpvalue>=significantline1)
+										HX1=(Cpvalue[p_amp.index]+r+H*(i-1)+cir.band*(i-1))*sin(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+										HY1=(Cpvalue[p_amp.index]+r+H*(i-1)+cir.band*(i-1))*cos(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+									}else{
+										if(LOG10){
+											significantline0=H*(-log10(threshold[ll-1])-Min)/(Max-Min)
+											significantline1=H*(-log10(threshold[ll])-Min)/(Max-Min)
+										}else{
+											significantline0=H*(threshold[ll-1]-Min)/(Max-Min)
+											significantline1=H*(threshold[ll]-Min)/(Max-Min)
+										}
+										p_amp.index <- which(Cpvalue>=significantline1 & Cpvalue < significantline0)
+										HX1=(Cpvalue[p_amp.index]+r+H*(i-1)+cir.band*(i-1))*sin(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+										HY1=(Cpvalue[p_amp.index]+r+H*(i-1)+cir.band*(i-1))*cos(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+									}
+								
+									if(is.null(signal.col)){
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[1],col=rep(rep(colx,N[i]),add[[i]])[p_amp.index])
+									}else{
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[1],col=signal.col[ll])
+									}
+								}
+						}
+					}
+				}
+				if(cir.chr==TRUE){
+					ticks1=1.07*(RR+cir.chr.h)*sin(2*pi*(ticks-round(band/2))/TotalN)
+					ticks2=1.07*(RR+cir.chr.h)*cos(2*pi*(ticks-round(band/2))/TotalN)
+					if(is.null(chr.labels)){
+						for(i in 1:length(ticks)){
+							angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+							text(ticks1[i],ticks2[i],chr.ori[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}else{
+						for(i in 1:length(ticks)){
+							angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+							text(ticks1[i],ticks2[i],chr.labels[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}
+				}else{
+					ticks1=(0.9*r)*sin(2*pi*(ticks-round(band/2))/TotalN)
+					ticks2=(0.9*r)*cos(2*pi*(ticks-round(band/2))/TotalN)
+					if(is.null(chr.labels)){
+						for(i in 1:length(ticks)){
+						angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+						text(ticks1[i],ticks2[i],chr.ori[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}else{
+						for(i in 1:length(ticks)){
+							angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+							text(ticks1[i],ticks2[i],chr.labels[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}
+				}
+			}
+			if(outward==FALSE){
+				if(cir.chr==TRUE){
+					# XLine=(2*cir.band+RR+cir.chr.h)*sin(2*pi*(1:TotalN)/TotalN)
+					# YLine=(2*cir.band+RR+cir.chr.h)*cos(2*pi*(1:TotalN)/TotalN)
+					# lines(XLine,YLine,lwd=1.5)
+
+					polygon.num <- 1000
+					for(k in 1:length(chr)){
+						if(k==1){
+							polygon.index <- seq(round(band/2)+1,-round(band/2)+max(pvalue.posN.list[[1]]), length=polygon.num)
+							X1chr=(2*cir.band+RR)*sin(2*pi*(polygon.index)/TotalN)
+							Y1chr=(2*cir.band+RR)*cos(2*pi*(polygon.index)/TotalN)
+							X2chr=(2*cir.band+RR+cir.chr.h)*sin(2*pi*(polygon.index)/TotalN)
+							Y2chr=(2*cir.band+RR+cir.chr.h)*cos(2*pi*(polygon.index)/TotalN)
+								if(is.null(chr.den.col)){
+									polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=rep(colx,ceiling(length(chr)/length(colx)))[k],border=rep(colx,ceiling(length(chr)/length(colx)))[k])	
+								}else{
+									if(cir.density){
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
+									}else{
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
+									}
+								}
+						}else{
+							polygon.index <- seq(1+round(band/2)+max(pvalue.posN.list[[k-1]]),-round(band/2)+max(pvalue.posN.list[[k]]), length=polygon.num)
+							X1chr=(2*cir.band+RR)*sin(2*pi*(polygon.index)/TotalN)
+							Y1chr=(2*cir.band+RR)*cos(2*pi*(polygon.index)/TotalN)
+							X2chr=(2*cir.band+RR+cir.chr.h)*sin(2*pi*(polygon.index)/TotalN)
+							Y2chr=(2*cir.band+RR+cir.chr.h)*cos(2*pi*(polygon.index)/TotalN)
+							if(is.null(chr.den.col)){
+								polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=rep(colx,ceiling(length(chr)/length(colx)))[k],border=rep(colx,ceiling(length(chr)/length(colx)))[k])	
+							}else{
+									if(cir.density){
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
+									}else{
+										polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
+									}
+							}	
+						}
+					}
+					if(cir.density){
+
+						segments(
+							(2*cir.band+RR)*sin(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							(2*cir.band+RR)*cos(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							(2*cir.band+RR+cir.chr.h)*sin(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							(2*cir.band+RR+cir.chr.h)*cos(2*pi*(pvalue.posN-round(band/2))/TotalN),
+							col=density.list$den.col, lwd=0.1
+						)
+						legend(
+							x=RR+4*cir.chr.h,
+							y=(RR+4*cir.chr.h)/2,
+							title="", legend=density.list$legend.y, pch=15, pt.cex = 3, col=density.list$legend.col,
+							cex=1, bty="n",
+							y.intersp=1,
+							x.intersp=1,
+							yjust=0.5, xjust=0, xpd=TRUE
+						)
+						
+					}
+					
+					if(cir.density){
+						circle.plot(myr=2*cir.band+RR+cir.chr.h,lwd=1.5,add=TRUE,col='grey')
+						circle.plot(myr=2*cir.band+RR,lwd=1.5,add=TRUE,col='grey')
+					}else{
+						circle.plot(myr=2*cir.band+RR+cir.chr.h,lwd=1.5,add=TRUE)
+						circle.plot(myr=2*cir.band+RR,lwd=1.5,add=TRUE)
+					}
+
+				}
+				
+				X=(-Cpvalue+r+H*i+cir.band*(i-1))*sin(2*pi*(pvalue.posN-round(band/2))/TotalN)
+				Y=(-Cpvalue+r+H*i+cir.band*(i-1))*cos(2*pi*(pvalue.posN-round(band/2))/TotalN)
+				points(X,Y,pch=19,cex=cex[1],col=rep(rep(colx,N[i]),add[[i]]))
+				
+								if(cir.legend==TRUE){
+					
+					#try to get the number after radix point
+					if((Max-Min)<=1) {
+						round.n=nchar(as.character(10^(-ceiling(-log10(Max)))))-1
+					}else{
+						round.n=2
+					}
+					segments(0,r+H*(i-1)+cir.band*(i-1),0,r+H*i+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					segments(0,r+H*(i-1)+cir.band*(i-1),H/20,r+H*(i-1)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-1)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0.75)+cir.band*(i-1),H/20,r+H*(i-0.75)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0.75)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0.5)+cir.band*(i-1),H/20,r+H*(i-0.5)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0.5)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0.25)+cir.band*(i-1),H/20,r+H*(i-0.25)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0.25)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					segments(0,r+H*(i-0)+cir.band*(i-1),H/20,r+H*(i-0)+cir.band*(i-1),col=cir.legend.col,lwd=1.5)
+					circle.plot(myr=r+H*(i-0)+cir.band*(i-1),lwd=0.5,add=TRUE,col='grey')
+					text(-r/15,r+H*(i-0.06)+cir.band*(i-1),round(Min+(Max-Min)*0,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.25)+cir.band*(i-1),round(Min+(Max-Min)*0.25,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.5)+cir.band*(i-1),round(Min+(Max-Min)*0.5,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.75)+cir.band*(i-1),round(Min+(Max-Min)*0.75,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+					text(-r/15,r+H*(i-0.94)+cir.band*(i-1),round(Min+(Max-Min)*1,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
+				}
+				
+				if(!is.null(threshold)){
+					if(sum(threshold!=0)==length(threshold)){
+					
+						for(thr in 1:length(threshold)){
+							significantline1=ifelse(LOG10, H*(-log10(threshold[thr])-Min)/(Max-Min), H*(threshold[thr]-Min)/(Max-Min))
+							#s1X=(significantline1+r+H*(i-1)+cir.band*(i-1))*sin(2*pi*(0:TotalN)/TotalN)
+							#s1Y=(significantline1+r+H*(i-1)+cir.band*(i-1))*cos(2*pi*(0:TotalN)/TotalN)
+							if(significantline1<H){
+								#lines(s1X,s1Y,type="l",col=threshold.col,lwd=threshold.col,lty=threshold.lty)
+								circle.plot(myr=(-significantline1+r+H*i+cir.band*(i-1)),col=threshold.col[thr],lwd=threshold.lwd[thr],lty=threshold.lty[thr])
+							}else{
+								warning(paste("No significant points for ",taxa[i]," pass the threshold level using threshold=",threshold[thr],"!",sep=""))
+							}
+						}
+						if(amplify==TRUE){
+							if(LOG10){
+								threshold <- sort(threshold)
+								significantline1=H*(-log10(max(threshold))-Min)/(Max-Min)
+							}else{
+								threshold <- sort(threshold, decreasing=TRUE)
+								significantline1=H*(min(threshold)-Min)/(Max-Min)
+							}
+							p_amp.index <- which(Cpvalue>=significantline1)
+							HX1=(-Cpvalue[p_amp.index]+r+H*i+cir.band*(i-1))*sin(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+							HY1=(-Cpvalue[p_amp.index]+r+H*i+cir.band*(i-1))*cos(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+							
+							#cover the points that exceed the threshold with the color "white"
+							points(HX1,HY1,pch=19,cex=cex[1],col="white")
+							
+								for(ll in 1:length(threshold)){
+									if(ll == 1){
+										if(LOG10){
+											significantline1=H*(-log10(threshold[ll])-Min)/(Max-Min)
+										}else{
+											significantline1=H*(threshold[ll]-Min)/(Max-Min)
+										}
+										p_amp.index <- which(Cpvalue>=significantline1)
+										HX1=(-Cpvalue[p_amp.index]+r+H*i+cir.band*(i-1))*sin(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+										HY1=(-Cpvalue[p_amp.index]+r+H*i+cir.band*(i-1))*cos(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+									}else{
+										if(LOG10){
+											significantline0=H*(-log10(threshold[ll-1])-Min)/(Max-Min)
+											significantline1=H*(-log10(threshold[ll])-Min)/(Max-Min)
+										}else{
+											significantline0=H*(threshold[ll-1]-Min)/(Max-Min)
+											significantline1=H*(threshold[ll]-Min)/(Max-Min)
+										}
+										p_amp.index <- which(Cpvalue>=significantline1 & Cpvalue < significantline0)
+										HX1=(-Cpvalue[p_amp.index]+r+H*i+cir.band*(i-1))*sin(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+										HY1=(-Cpvalue[p_amp.index]+r+H*i+cir.band*(i-1))*cos(2*pi*(pvalue.posN[p_amp.index]-round(band/2))/TotalN)
+									
+									}
+								
+									if(is.null(signal.col)){
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[1],col=rep(rep(colx,N[i]),add[[i]])[p_amp.index])
+									}else{
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[1],col=signal.col[ll])
+									}
+								}
+						}
+					}
+				}
+				
+				if(cir.chr==TRUE){
+					ticks1=1.1*(2*cir.band+RR)*sin(2*pi*(ticks-round(band/2))/TotalN)
+					ticks2=1.1*(2*cir.band+RR)*cos(2*pi*(ticks-round(band/2))/TotalN)
+					if(is.null(chr.labels)){
+						for(i in 1:length(ticks)){
+						  angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+						  text(ticks1[i],ticks2[i],chr.ori[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}else{
+						for(i in 1:length(ticks)){
+							angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+							text(ticks1[i],ticks2[i],chr.labels[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}
+				}else{
+					ticks1=1.0*(RR+cir.band)*sin(2*pi*(ticks-round(band/2))/TotalN)
+					ticks2=1.0*(RR+cir.band)*cos(2*pi*(ticks-round(band/2))/TotalN)
+					if(is.null(chr.labels)){
+						for(i in 1:length(ticks)){
+						
+							#adjust the angle of labels of circle plot
+							angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+							text(ticks1[i],ticks2[i],chr.ori[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}else{
+						for(i in 1:length(ticks)){
+							angle=360*(1-(ticks-round(band/2))[i]/TotalN)
+							text(ticks1[i],ticks2[i],chr.labels[i],srt=angle,font=2,cex=cex.axis)
+						}
+					}	
+				}
+			}
+		}
+		if(file.output) dev.off()
+		#print("Circular-Manhattan has been finished!",quote=F)
+	}
+
+	if("m" %in% plot.type){
+		if(multracks==FALSE){
+			#print("Starting Rectangular-Manhattan plot!",quote=F)
+			for(i in 1:R){
+				colx=col[i,]
+				colx=colx[!is.na(colx)]
+				if(verbose)	print(paste("Rectangular_Manhattan Plotting ",taxa[i],"...",sep=""))
+					if(file.output){
+						if(file=="jpg")	jpeg(paste("Manhattan.",taxa[i],".jpg",sep=""), width = 14*dpi,height=5*dpi,res=dpi,quality = 100)
+						if(file=="pdf")	pdf(paste("Manhattan.",taxa[i],".pdf",sep=""), width = 15,height=6)
+						if(file=="tiff")	tiff(paste("Manhattan.",taxa[i],".tiff",sep=""), width = 14*dpi,height=5*dpi,res=dpi)
+						par(mar = c(5,6,4,3),xaxs=xaxs,yaxs=yaxs,xpd=TRUE)
+					}
+					if(!file.output){
+						if(is.null(dev.list()))	dev.new(width = 15, height = 6)
+						par(xpd=TRUE)
+					}
+					
+					pvalue=pvalueT[,i]
+					logpvalue=logpvalueT[,i]
+					if(is.null(ylim)){
+						if(!is.null(threshold)){
+							if(sum(threshold!=0)==length(threshold)){
+								if(LOG10 == TRUE){
+									Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))),ceiling(-log10(min(threshold))))
+									Min <- 0
+								}else{
+									Max=max(ceiling(max(pvalue[pvalue!=Inf])),max(threshold))
+									if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]),max(threshold))
+									Min <- min(floor(min(pvalue[pvalue!=Inf])),min(threshold))
+									if(abs(Min)<=1)	Min=min(min(pvalue[pvalue!=Inf]),min(threshold))
+								}
+							}else{
+								if(LOG10){
+									Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))))
+									Min<-0
+								}else{
+									Max=ceiling(max(pvalue[pvalue!=Inf]))
+									if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]))
+									Min<-floor(min(pvalue[pvalue!=Inf]))
+									if(abs(Min)<=1)	Min=min(pvalue[pvalue!=Inf])
+									# }else{
+										# Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+									# }
+								}
+							}
+						}else{
+							if(LOG10){
+								Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))))
+								Min=0
+							}else{
+								Max=ceiling(max(pvalue[pvalue!=Inf]))
+								if(abs(Max)<=1)	Max=max(pvalue[pvalue!=Inf])
+								Min=floor(min(pvalue[pvalue!=Inf]))
+								if(abs(Min)<=1)	Min=min(pvalue[pvalue!=Inf])
+								# }else{
+									# Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+								# }
+							}
+						}
+						if((Max-Min)<=1){
+							if(cir.density){
+								plot(pvalue.posN,logpvalue,pch=pch,cex=cex[2],col=rep(rep(colx,N[i]),add[[i]]),xlim=c(0,1.01*max(pvalue.posN)),ylim=c(Min-(Max-Min)/den.fold, Max+10^(-ceiling(-log10(abs(Max))))),ylab=ylab,
+									cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main=taxa)
+							}else{
+								plot(pvalue.posN,logpvalue,pch=pch,cex=cex[2],col=rep(rep(colx,N[i]),add[[i]]),xlim=c(0,max(pvalue.posN)),ylim=c(Min,Max+10^(-ceiling(-log10(abs(Max))))),ylab=ylab,
+								cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main=taxa)
+							}
+						}else{
+							if(cir.density){
+								plot(pvalue.posN,logpvalue,pch=pch,cex=cex[2],col=rep(rep(colx,N[i]),add[[i]]),xlim=c(0,1.01*max(pvalue.posN)),ylim=c(Min-(Max-Min)/den.fold,Max),ylab=ylab,
+								cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main=taxa)
+							}else{
+								plot(pvalue.posN,logpvalue,pch=pch,cex=cex[2],col=rep(rep(colx,N[i]),add[[i]]),xlim=c(0,max(pvalue.posN)),ylim=c(Min,Max),ylab=ylab,
+								cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main=taxa)
+							}
+						}
+					}else{
+						Max <- max(ylim)
+						Min <- min(ylim)
+						if(cir.density){
+							plot(pvalue.posN[logpvalue>=min(ylim)],logpvalue[logpvalue>=min(ylim)],pch=pch,cex=cex[2],col=rep(rep(colx,N[i]),add[[i]])[logpvalue>=min(ylim)],xlim=c(0,1.01*max(pvalue.posN)),ylim=c(min(ylim)-(Max-Min)/den.fold, max(ylim)),ylab=ylab,
+							cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main=taxa)
+						}else{
+							plot(pvalue.posN[logpvalue>=min(ylim)],logpvalue[logpvalue>=min(ylim)],pch=pch,cex=cex[2],col=rep(rep(colx,N[i]),add[[i]])[logpvalue>=min(ylim)],xlim=c(0,max(pvalue.posN)),ylim=ylim,ylab=ylab,
+							cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main=taxa)
+						}
+					}
+					Max1 <- Max
+					Min1 <- Min
+					if(abs(Max) <= 1) Max <- round(Max, ceiling(-log10(abs(Max))))
+					if(abs(Min) <= 1) Min <- round(Min, ceiling(-log10(abs(Min))))
+					if(is.null(chr.labels)){
+						axis(1, at=c(0,ticks),cex.axis=cex.axis,font=2,labels=c("Chr",chr.ori))
+					}else{
+						axis(1, at=c(0,ticks),cex.axis=cex.axis,font=2,labels=c("Chr",chr.labels))
+					}
+					if(is.null(ylim)){
+						if((Max-Min)>1){
+							axis(2,at=seq(Min,(Max),ceiling((Max-Min)/10)),cex.axis=cex.axis,font=2,labels=round(seq(Min,(Max),ceiling((Max-Min)/10)), 2))
+							legend.y <- tail(round(seq(Min,(Max),ceiling((Max-Min)/10)), 2), 1)
+						}else{
+							axis(2,at=seq(Min,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))),cex.axis=cex.axis,font=2,labels=seq(Min,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))))
+							legend.y <- tail(seq(Min,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))), 1)
+						}
+					}else{
+						if(ylim[2]>1){
+							axis(2,at=seq(min(ylim),ylim[2],ceiling((ylim[2])/10)),cex.axis=cex.axis,font=2,labels=round(seq(min(ylim),(ylim[2]),ceiling((ylim[2])/10)), 2))
+							legend.y <- tail(ylim[2], 1)
+						}else{
+							axis(2,at=seq(min(ylim),ylim[2],10^(-ceiling(-log10(max(abs(ylim)))))),cex.axis=cex.axis,font=2,labels=seq(min(ylim),ylim[2],10^(-ceiling(-log10(max(abs(ylim)))))))
+							legend.y <- tail(ylim[2], 1)
+						}
+					}
+					if(!is.null(threshold)){
+						if(sum(threshold!=0)==length(threshold)){
+							for(thr in 1:length(threshold)){
+								h <- ifelse(LOG10, -log10(threshold[thr]), threshold[thr])
+								# print(h)
+								# print(threshold.col[thr])
+								# print(threshold.lty[thr])
+								# print(threshold.lwd[thr])
+								par(xpd=FALSE); abline(h=h,col=threshold.col[thr],lty=threshold.lty[thr],lwd=threshold.lwd[thr]); par(xpd=TRUE)
+							}
+							if(amplify == TRUE){
+								if(LOG10){
+									threshold <- sort(threshold)
+									sgline1=-log10(max(threshold))
+								}else{
+									threshold <- sort(threshold, decreasing=TRUE)
+									sgline1=min(threshold)
+								}
+
+								sgindex=which(logpvalue>=sgline1)
+								HY1=logpvalue[sgindex]
+								HX1=pvalue.posN[sgindex]
+								
+								#cover the points that exceed the threshold with the color "white"
+								points(HX1,HY1,pch=pch,cex=cex[2],col="white")
+								
+								for(ll in 1:length(threshold)){
+									if(ll == 1){
+										if(LOG10){
+											sgline1=-log10(threshold[ll])
+										}else{
+											sgline1=threshold[ll]
+										}
+										sgindex=which(logpvalue>=sgline1)
+										HY1=logpvalue[sgindex]
+										HX1=pvalue.posN[sgindex]
+									}else{
+										if(LOG10){
+											sgline0=-log10(threshold[ll-1])
+											sgline1=-log10(threshold[ll])
+										}else{
+											sgline0=threshold[ll-1]
+											sgline1=threshold[ll]
+										}
+										sgindex=which(logpvalue>=sgline1 & logpvalue < sgline0)
+										HY1=logpvalue[sgindex]
+										HX1=pvalue.posN[sgindex]
+									}
+
+									if(is.null(signal.col)){
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[2],col=rep(rep(colx,N[i]),add[[i]])[sgindex])
+									}else{
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[2],col=signal.col[ll])
+									}
+									
+								}
+							}
+						}
+					}
+					if(is.null(ylim)){ymin <- Min1}else{ymin <- min(ylim)}
+					if(cir.density){
+						for(yll in 1:length(pvalue.posN.list)){
+							polygon(c(min(pvalue.posN.list[[yll]]), min(pvalue.posN.list[[yll]]), max(pvalue.posN.list[[yll]]), max(pvalue.posN.list[[yll]])), 
+								c(ymin-0.5*(Max-Min)/den.fold, ymin-1.5*(Max-Min)/den.fold, 
+								ymin-1.5*(Max-Min)/den.fold, ymin-0.5*(Max-Min)/den.fold), 
+								col="grey", border="grey")
+						}
+						
+						segments(
+							pvalue.posN,
+							ymin-0.5*(Max-Min)/den.fold,
+							pvalue.posN,
+							ymin-1.5*(Max-Min)/den.fold,
+							col=density.list$den.col, lwd=0.1
+						)
+						legend(
+							x=max(pvalue.posN)+band,
+							y=legend.y,
+							title="", legend=density.list$legend.y, pch=15, pt.cex = 2.5, col=density.list$legend.col,
+							cex=0.8, bty="n",
+							y.intersp=1,
+							x.intersp=1,
+							yjust=1, xjust=0, xpd=TRUE
+						)
+						
+					}
+				if(box) box()
+				#if(!is.null(threshold) & (length(grep("FarmCPU",taxa[i])) != 0))	abline(v=which(pvalueT[,i] < min(threshold)/max(dim(Pmap))),col="grey",lty=2,lwd=signal.line)
+				if(file.output)  dev.off()
+			}
+			#print("Rectangular-Manhattan has been finished!",quote=F)
+		}else{
+			#print("Starting Rectangular-Manhattan plot!",quote=F)
+			#print("Plotting in multiple tracks!",quote=F)
+			if(file.output){
+				if(file=="jpg")	jpeg(paste("Multracks.Rectangular-Manhattan.",paste(taxa,collapse="."),".jpg",sep=""), width = 14*dpi,height=5*dpi*R,res=dpi,quality = 100)
+				if(file=="pdf")	pdf(paste("Multracks.Rectangular-Manhattan.",paste(taxa,collapse="."),".pdf",sep=""), width = 15,height=6*R)
+				if(file=="tiff")	tiff(paste("Multracks.Rectangular-Manhattan.",paste(taxa,collapse="."),".tiff",sep=""), width = 14*dpi,height=5*dpi*R,res=dpi)
+				par(mfcol=c(R,1),mar=c(0, 6+(R-1)*2, 0, 2),oma=c(4,0,4,0),xaxs=xaxs,yaxs=yaxs,xpd=TRUE)
+			}
+			if(!file.output){
+				if(is.null(dev.list()))	dev.new(width = 15, height = 6)
+				par(xpd=TRUE)
+			}
+			for(i in 1:R){
+				if(verbose)	print(paste("Multracks_Rectangular Plotting ",taxa[i],"...",sep=""))
+				colx=col[i,]
+				colx=colx[!is.na(colx)]
+				pvalue=pvalueT[,i]
+				logpvalue=logpvalueT[,i]
+				if(is.null(ylim)){
+					if(!is.null(threshold)){
+						if(sum(threshold!=0)==length(threshold)){
+							if(LOG10){
+								Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))),-log10(min(threshold)))
+								Min <- 0
+							}else{
+								Max=max(ceiling(max(pvalue[pvalue!=Inf])),max(threshold))
+								if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]),max(threshold))
+								Min<-min(floor(min(pvalue[pvalue!=Inf])),min(threshold))
+								if(abs(Min)<=1)	Min=min(min(pvalue[pvalue!=Inf]),min(threshold))
+							}
+						}else{
+							if(LOG10){
+								Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))))
+								Min<-0
+							}else{
+								Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+								if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]))
+								Min=min(floor(min(pvalue[pvalue!=Inf])))
+								if(abs(Min)<=1)	Min=min(min(pvalue[pvalue!=Inf]))
+								# }else{
+									# Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+								# }
+							}	
+						}
+					}else{
+						if(LOG10){
+							Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))))
+							Min=0
+						}else{
+							Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+							if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]))
+							Min <- min(ceiling(min(pvalue[pvalue!=Inf])))
+							if(abs(Min)<=1)	Min=min(min(pvalue[pvalue!=Inf]))
+							# }else{
+								# Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+							# }
+						}
+					}
+					xn <- ifelse(R == 1, R, R * 2/3)
+					if((Max-Min)<=1){
+						plot(pvalue.posN,logpvalue,pch=pch,cex=cex[2]*xn,col=rep(rep(colx,N[i]),add[[i]]),xlim=c(0,max(pvalue.posN)+band),ylim=c(Min,Max+10^(-ceiling(-log10(abs(Max))))),ylab=ylab,
+							cex.axis=cex.axis*xn,cex.lab=2*xn,font=2,axes=FALSE)
+					}else{
+						plot(pvalue.posN,logpvalue,pch=pch,cex=cex[2]*xn,col=rep(rep(colx,N[i]),add[[i]]),xlim=c(0,max(pvalue.posN)+band),ylim=c(Min,Max),ylab=ylab,
+							cex.axis=cex.axis*xn,cex.lab=2*xn,font=2,axes=FALSE)
+					}
+				}else{
+					xn <- ifelse(R == 1, R, R * 2/3)
+					Max <- max(ylim)
+					Min <- min(ylim)
+					plot(pvalue.posN[logpvalue>=min(ylim)],logpvalue[logpvalue>=min(ylim)],pch=pch,cex=cex[2]*xn,col=rep(rep(colx,N[i]),add[[i]])[logpvalue>=min(ylim)],xlim=c(0,max(pvalue.posN)+band),ylim=ylim,ylab=ylab,
+						cex.axis=cex.axis*xn,cex.lab=2*xn,font=2,axes=FALSE)
+				}
+				Max1 <- Max
+				Min1 <- Min
+				if(abs(Max) <= 1) Max <- round(Max, ceiling(-log10(abs(Max))))
+				if(abs(Min) <= 1) Min <- round(Min, ceiling(-log10(abs(Min))))
+				
+				#add the names of traits on plot 
+				if((Max-Min)<=1){
+					text(ticks[1],Max+10^(-ceiling(-log10(abs(Max)))),labels=taxa[i],adj=0,font=3,cex=xn)
+				}else{
+					text(ticks[1],Max+1,labels=taxa[i],adj=0,font=3,cex=xn)
+				}
+				if(i == R){
+					if(is.null(chr.labels)){
+						axis(1, at=c(0,ticks),cex.axis=cex.axis*xn,font=2,labels=c("Chr",chr.ori),padj=(xn-1)/2)
+					}else{
+						axis(1, at=c(0,ticks),cex.axis=cex.axis*xn,font=2,labels=c("Chr",chr.labels),padj=(xn-1)/2)
+					}
+				}
+				#if(i==1) mtext("Manhattan plot",side=3,padj=-1,font=2,cex=xn)
+				if(is.null(ylim)){
+					if((Max-Min)>1){
+						axis(2,at=seq(Min,(Max),ceiling((Max-Min)/10)),cex.axis=cex.axis*xn,font=2,labels=round(seq(Min,(Max),ceiling((Max-Min)/10)), 2))
+					}else{
+						axis(2,at=seq(Min,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))),cex.axis=cex.axis*xn,font=2,labels=seq(0,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))))
+					}
+				}else{
+					if(ylim[2]>1){
+						axis(2,at=seq(min(ylim),(ylim[2]),ceiling((ylim[2])/10)),cex.axis=cex.axis*xn,font=2,labels=round(seq(min(ylim),(ylim[2]),ceiling((ylim[2])/10)), 2))
+					}else{
+						axis(2,at=seq(min(ylim),ylim[2],10^(-ceiling(-log10(max(abs(ylim)))))),cex.axis=cex.axis*xn,font=2,labels=seq(min(ylim),ylim[2],10^(-ceiling(-log10(max(abs(ylim)))))))
+					}
+				}
+				if(!is.null(threshold)){
+					if(sum(threshold!=0)==length(threshold)){
+						for(thr in 1:length(threshold)){
+							h <- ifelse(LOG10, -log10(threshold[thr]), threshold[thr])
+							par(xpd=FALSE); abline(h=h,col=threshold.col[thr],lwd=threshold.lwd[thr],lty=threshold.lty[thr]); par(xpd=TRUE)
+						}
+						if(amplify==TRUE){
+								if(LOG10){
+									threshold <- sort(threshold)
+									sgline1=-log10(max(threshold))
+								}else{
+									threshold <- sort(threshold, decreasing=TRUE)
+									sgline1=min(threshold)
+								}
+								sgindex=which(logpvalue>=sgline1)
+								HY1=logpvalue[sgindex]
+								HX1=pvalue.posN[sgindex]
+								
+								#cover the points that exceed the threshold with the color "white"
+								points(HX1,HY1,pch=pch,cex=cex[2]*xn,col="white")
+								
+								for(ll in 1:length(threshold)){
+									if(ll == 1){
+										if(LOG10){
+											sgline1=-log10(threshold[ll])
+										}else{
+											sgline1=threshold[ll]
+										}
+										sgindex=which(logpvalue>=sgline1)
+										HY1=logpvalue[sgindex]
+										HX1=pvalue.posN[sgindex]
+									}else{
+										if(LOG10){
+											sgline0=-log10(threshold[ll-1])
+											sgline1=-log10(threshold[ll])
+										}else{
+											sgline0=threshold[ll-1]
+											sgline1=threshold[ll]
+										}
+										sgindex=which(logpvalue>=sgline1 & logpvalue < sgline0)
+										HY1=logpvalue[sgindex]
+										HX1=pvalue.posN[sgindex]
+									}
+
+									if(is.null(signal.col)){
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[2]*xn,col=rep(rep(colx,N[i]),add[[i]])[sgindex])
+									}else{
+										points(HX1,HY1,pch=signal.pch[ll],cex=signal.cex[ll]*cex[2]*xn,col=signal.col[ll])
+									}
+									
+								}
+						}
+					}
+				}
+				#if(!is.null(threshold) & (length(grep("FarmCPU",taxa[i])) != 0))	abline(v=which(pvalueT[,i] < min(threshold)/max(dim(Pmap))),col="grey",lty=2,lwd=signal.line)
+			}
+			
+			#add the labels of X-axis
+			#mtext(xlab,side=1,padj=2.5,font=2,cex=R*2/3)
+			if(file.output) dev.off()
+			
+			if(file.output){
+				if(file=="jpg")	jpeg(paste("Multraits.Rectangular-Manhattan.",paste(taxa,collapse="."),".jpg",sep=""), width = 14*dpi,height=5*dpi,res=dpi,quality = 100)
+				if(file=="pdf")	pdf(paste("Multraits.Rectangular-Manhattan.",paste(taxa,collapse="."),".pdf",sep=""), width = 15,height=6)
+				if(file=="tiff")	tiff(paste("Multraits.Rectangular-Manhattan.",paste(taxa,collapse="."),".tiff",sep=""), width = 14*dpi,height=5*dpi,res=dpi)
+				par(mar = c(5,6,4,3),xaxs=xaxs,yaxs=yaxs,xpd=TRUE)
+			}
+			if(!file.output){
+				if(is.null(dev.list()))	dev.new(width = 15, height = 6)
+				par(xpd=TRUE)
+			}
+			
+			pvalue <- as.vector(Pmap[,3:(R+2)])
+			if(is.null(ylim)){
+				if(!is.null(threshold)){
+					if(sum(threshold!=0)==length(threshold)){
+						if(LOG10){
+							Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))),-log10(min(threshold)))
+							Min<-0
+						}else{
+							Max=max(ceiling(max(pvalue[pvalue!=Inf])),max(threshold))
+							if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]),max(threshold))
+							Min <- min(floor(min(pvalue[pvalue!=Inf])),min(threshold))
+							if(abs(Min)<=1)	Min=min(min(pvalue[pvalue!=Inf]),min(threshold))
+						}
+					}else{
+						if(LOG10){
+							Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))))
+							Min=0
+						}else{
+							Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+							if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]))
+							Min<- min(floor(min(pvalue[pvalue!=Inf])))
+							if(abs(Min)<=1)	Min=min(min(pvalue[pvalue!=Inf]))
+							# }else{
+								# Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+							# }
+						}	
+					}
+				}else{
+					if(LOG10){
+						Max=max(ceiling(-log10(min(pvalue[pvalue!=0]))))
+						Min=0
+					}else{
+						Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+						
+						#{
+						if(abs(Max)<=1)	Max=max(max(pvalue[pvalue!=Inf]))
+						Min<- min(floor(min(pvalue[pvalue!=Inf])))
+						if(abs(Min)<=1)	Min=min(min(pvalue[pvalue!=Inf]))
+						
+						# }else{
+							# Max=max(ceiling(max(pvalue[pvalue!=Inf])))
+						# }
+					}
+				}
+				if((Max-Min)<=1){
+					if(cir.density){
+						plot(NULL,xlim=c(0,1.01*max(pvalue.posN)),ylim=c(Min-(Max-Min)/den.fold, Max+10^(-ceiling(-log10(abs(Max))))),ylab=ylab,
+							cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main="")
+					}else{
+						plot(NULL,xlim=c(0,max(pvalue.posN)),ylim=c(Min,Max+10^(-ceiling(-log10(abs(Max))))),ylab=ylab,
+							cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main="")
+					}
+				}else{
+					if(cir.density){
+						plot(NULL,xlim=c(0,1.01*max(pvalue.posN)),ylim=c(Min-(Max-Min)/den.fold,Max),ylab=ylab,
+							cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main="")
+					}else{
+						plot(NULL,xlim=c(0,max(pvalue.posN)),ylim=c(Min,Max),ylab=ylab,
+							cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main="")
+					}
+				}
+			}else{
+				Max <- max(ylim)
+				Min <- min(ylim)
+				if(cir.density){
+					plot(NULL,xlim=c(0,1.01*max(pvalue.posN)),ylim=c(min(ylim)-Max/den.fold,Max),ylab=ylab,
+						cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main="Manhattan plot of")
+				}else{
+					plot(NULL,xlim=c(0,max(pvalue.posN)),ylim=ylim,ylab=ylab,
+						cex.axis=cex.axis,cex.lab=2,font=2,axes=FALSE,xlab=xlab,main="Manhattan plot of")
+				}
+			}
+			Max1 <- Max
+			Min1 <- Min
+			if(abs(Max) <= 1) Max <- round(Max, ceiling(-log10(abs(Max))))
+			if(abs(Min) <= 1) Min <- round(Min, ceiling(-log10(abs(Min))))
+			legend("topleft",taxa,col=t(col)[1:R],pch=19,text.font=6,box.col=NA)
+			if(is.null(chr.labels)){
+				axis(1, at=c(0,ticks),cex.axis=cex.axis,font=2,labels=c("Chr",chr.ori))
+			}else{
+				axis(1, at=c(0,ticks),cex.axis=cex.axis,font=2,labels=c("Chr",chr.labels))
+			}
+			if(is.null(ylim)){
+				if((Max-Min)>1){
+					#print(seq(0,(Max+1),ceiling((Max+1)/10)))
+					axis(2,at=seq(Min,(Max),ceiling((Max-Min)/10)),cex.axis=cex.axis,font=2,labels=round(seq(Min,(Max),ceiling((Max-Min)/10)),2))
+					legend.y <- tail(round(seq(0,(Max),ceiling((Max)/10)),2), 1)
+				}else{
+					axis(2,at=seq(Min,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))),cex.axis=cex.axis,font=2,labels=seq(0,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))))
+					legend.y <- tail(seq(0,Max+10^(-ceiling(-log10(abs(Max)))),10^(-ceiling(-log10(max(abs(c(Max, Min))))))), 1)
+				}
+			}else{
+				if(ylim[2]>1){
+					axis(2,at=seq(min(ylim),(ylim[2]),ceiling((ylim[2])/10)),cex.axis=cex.axis,font=2,labels=round(seq(min(ylim),(ylim[2]),ceiling((ylim[2])/10)),2))
+					legend.y <- tail(ylim[2], 1)
+				}else{
+					axis(2,at=seq(min(ylim),ylim[2],10^(-ceiling(-log10(max(abs(ylim)))))),cex.axis=cex.axis,font=2,labels=seq(min(ylim),ylim[2],10^(-ceiling(-log10(max(abs(ylim)))))))
+					legend.y <- tail(ylim[2], 1)
+				}
+			}
+			do <- TRUE
+			sam.index <- list()
+			for(l in 1:R){
+				sam.index[[l]] <- 1:nrow(Pmap)
+			}
+			
+			#change the sample number according to Pmap
+			sam.num <- ceiling(nrow(Pmap)/30)
+			if(verbose)	print("Multraits_Rectangular Plotting...")
+			while(do){
+				for(i in 1:R){
+					if(length(sam.index[[i]]) < sam.num){
+						plot.index <- sam.index[[i]]
+					}else{
+						plot.index <- sample(sam.index[[i]], sam.num, replace=FALSE)
+					}
+					sam.index[[i]] <- sam.index[[i]][-which(sam.index[[i]] %in% plot.index)]
+					logpvalue=logpvalueT[plot.index,i]
+					if(!is.null(ylim)){indexx <- logpvalue>=min(ylim)}else{indexx <- 1:length(logpvalue)}
+					points(pvalue.posN[plot.index][indexx],logpvalue[indexx],pch=pch,cex=cex[2],col=rgb(col2rgb(t(col)[i])[1], col2rgb(t(col)[i])[2], col2rgb(t(col)[i])[3], 100, maxColorValue=255))
+					#if(!is.null(threshold) & (length(grep("FarmCPU",taxa[i])) != 0))	abline(v=which(pvalueT[,i] < min(threshold)/max(dim(Pmap))),col="grey",lty=2,lwd=signal.line)
+				}
+				if(length(sam.index[[i]]) == 0) do <- FALSE
+			}
+			
+			# for(i in 1:R){
+				# logpvalue=logpvalueT[,i]
+				# points(pvalue.posN,logpvalue,pch=pch,cex=cex[2],col=t(col)[i])
+			# }
+			
+			if(!is.null(threshold)){
+				if(sum(threshold!=0)==length(threshold)){
+					for(thr in 1:length(threshold)){
+						h <- ifelse(LOG10, -log10(threshold[thr]), threshold[thr])
+						par(xpd=FALSE); abline(h=h,col=threshold.col[thr],lwd=threshold.lwd[thr],lty=threshold.lty[thr]); par(xpd=TRUE)
+					}
+				}
+			}
+			if(is.null(ylim)){ymin <- Min1}else{ymin <- min(ylim)}
+			if(cir.density){
+						for(yll in 1:length(pvalue.posN.list)){
+							polygon(c(min(pvalue.posN.list[[yll]]), min(pvalue.posN.list[[yll]]), max(pvalue.posN.list[[yll]]), max(pvalue.posN.list[[yll]])), 
+								c(ymin-0.5*(Max-Min)/den.fold, ymin-1.5*(Max-Min)/den.fold, 
+								ymin-1.5*(Max-Min)/den.fold, ymin-0.5*(Max-Min)/den.fold), 
+								col="grey", border="grey")
+						}
+						
+						segments(
+							pvalue.posN,
+							ymin-0.5*(Max-Min)/den.fold,
+							pvalue.posN,
+							ymin-1.5*(Max-Min)/den.fold,
+							col=density.list$den.col, lwd=0.1
+						)
+						legend(
+							x=max(pvalue.posN)+band,
+							y=legend.y,
+							title="", legend=density.list$legend.y, pch=15, pt.cex = 2.5, col=density.list$legend.col,
+							cex=0.8, bty="n",
+							y.intersp=1,
+							x.intersp=1,
+							yjust=1, xjust=0, xpd=TRUE
+						)
+						
+			}
+			if(file.output) dev.off()
+			
+		}
+	}
+		
+	if("q" %in% plot.type){
+		#print("Starting QQ-plot!",quote=F)
+		if(multracks){
+			if(file.output){
+				if(file=="jpg")	jpeg(paste("Multracks.QQplot.",paste(taxa,collapse="."),".jpg",sep=""), width = R*2.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
+				if(file=="pdf")	pdf(paste("Multracks.QQplot.",paste(taxa,collapse="."),".pdf",sep=""), width = R*2.5,height=5.5)
+				if(file=="tiff")	tiff(paste("Multracks.QQplot.",paste(taxa,collapse="."),".tiff",sep=""), width = R*2.5*dpi,height=5.5*dpi,res=dpi)
+				par(mfcol=c(1,R),mar = c(0,1,4,1.5),oma=c(3,5,0,0),xpd=TRUE)
+			}else{
+				if(is.null(dev.list()))	dev.new(width = 2.5*R, height = 5.5)
+				par(xpd=TRUE)
+			}
+			for(i in 1:R){
+				if(verbose)	print(paste("Multracks_QQ Plotting ",taxa[i],"...",sep=""))		
+				P.values=as.numeric(Pmap[,i+2])
+				P.values=P.values[!is.na(P.values)]
+				if(LOG10){
+					P.values=P.values[P.values>0]
+					P.values=P.values[P.values<=1]
+					N=length(P.values)
+					P.values=P.values[order(P.values)]
+				}else{
+					N=length(P.values)
+					P.values=P.values[order(P.values,decreasing=TRUE)]
+				}
+				p_value_quantiles=(1:length(P.values))/(length(P.values)+1)
+				log.Quantiles <- -log10(p_value_quantiles)
+				if(LOG10){
+					log.P.values <- -log10(P.values)
+				}else{
+					log.P.values <- P.values
+				}
+				
+				#calculate the confidence interval of QQ-plot
+				if(conf.int){
+					N1=length(log.Quantiles)
+					c95 <- rep(NA,N1)
+					c05 <- rep(NA,N1)
+					for(j in 1:N1){
+						xi=ceiling((10^-log.Quantiles[j])*N)
+						if(xi==0)xi=1
+						c95[j] <- qbeta(0.95,xi,N-xi+1)
+						c05[j] <- qbeta(0.05,xi,N-xi+1)
+					}
+					index=length(c95):1
+				}else{
+					c05 <- 1
+					c95 <- 1
+				}
+				
+				YlimMax <- max(floor(max(max(-log10(c05)), max(-log10(c95)))+1), floor(max(log.P.values)+1))
+				plot(NULL, xlim = c(0,floor(max(log.Quantiles)+1)), axes=FALSE, cex.axis=cex.axis, cex.lab=1.2,ylim=c(0,YlimMax),xlab ="", ylab="", main = taxa[i])
+				axis(1, at=seq(0,floor(max(log.Quantiles)+1),ceiling((max(log.Quantiles)+1)/10)), labels=seq(0,floor(max(log.Quantiles)+1),ceiling((max(log.Quantiles)+1)/10)), cex.axis=cex.axis)
+				axis(2, at=seq(0,YlimMax,ceiling(YlimMax/10)), labels=seq(0,YlimMax,ceiling(YlimMax/10)), cex.axis=cex.axis)
+				
+				#plot the confidence interval of QQ-plot
+				if(conf.int)	polygon(c(log.Quantiles[index],log.Quantiles),c(-log10(c05)[index],-log10(c95)),col=conf.int.col,border=conf.int.col)
+				
+				if(!is.null(threshold.col)){par(xpd=FALSE); abline(a = 0, b = 1, col = threshold.col[1],lwd=2); par(xpd=TRUE)}
+				points(log.Quantiles, log.P.values, col = col[1],pch=19,cex=cex[3])
+				if(!is.null(threshold)){
+					if(sum(threshold!=0)==length(threshold)){
+						thre.line=-log10(min(threshold))
+						if(amplify==TRUE){
+							thre.index=which(log.P.values>=thre.line)
+							if(length(thre.index)!=0){
+							
+								#cover the points that exceed the threshold with the color "white"
+								points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,cex=cex[3])
+								if(is.null(signal.col)){
+									points(log.Quantiles[thre.index],log.P.values[thre.index],col = col[1],pch=signal.pch[1],cex=signal.cex[1])
+								}else{
+									points(log.Quantiles[thre.index],log.P.values[thre.index],col = signal.col[1],pch=signal.pch[1],cex=signal.cex[1])
+								}
+							}
+						}
+					}
+				}
+			}
+			if(box)	box()
+			if(file.output) dev.off()
+			if(R > 1){
+				signal.col <- NULL
+				if(file.output){
+					if(file=="jpg")	jpeg(paste("Multraits.QQplot.",paste(taxa,collapse="."),".jpg",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
+					if(file=="pdf")	pdf(paste("Multraits.QQplot.",paste(taxa,collapse="."),".pdf",sep=""), width = 5.5,height=5.5)
+					if(file=="tiff")	tiff(paste("Multraits.QQplot.",paste(taxa,collapse="."),".tiff",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi)
+					par(mar = c(5,5,4,2),xpd=TRUE)
+				}else{
+					dev.new(width = 5.5, height = 5.5)
+					par(xpd=TRUE)
+				}
+				p_value_quantiles=(1:nrow(Pmap))/(nrow(Pmap)+1)
+				log.Quantiles <- -log10(p_value_quantiles)
+											
+				# calculate the confidence interval of QQ-plot
+				if((i == 1) & conf.int){
+					N1=length(log.Quantiles)
+					c95 <- rep(NA,N1)
+					c05 <- rep(NA,N1)
+					for(j in 1:N1){
+						xi=ceiling((10^-log.Quantiles[j])*N)
+						if(xi==0)xi=1
+						c95[j] <- qbeta(0.95,xi,N-xi+1)
+						c05[j] <- qbeta(0.05,xi,N-xi+1)
+					}
+					index=length(c95):1
+				}
+				
+				if(!conf.int){c05 <- 1; c95 <- 1}
+				
+				Pmap.min <- Pmap[,3:(R+2)]
+				YlimMax <- max(floor(max(max(-log10(c05)), max(-log10(c95)))+1), -log10(min(Pmap.min[Pmap.min > 0])))
+				plot(NULL, xlim = c(0,floor(max(log.Quantiles)+1)), axes=FALSE, cex.axis=cex.axis, cex.lab=1.2,ylim=c(0, floor(YlimMax+1)),xlab =expression(Expected~~-log[10](italic(p))), ylab = expression(Observed~~-log[10](italic(p))), main = "QQplot")
+				legend("topleft",taxa,col=t(col)[1:R],pch=19,text.font=6,box.col=NA)
+				axis(1, at=seq(0,floor(max(log.Quantiles)+1),ceiling((max(log.Quantiles)+1)/10)), labels=seq(0,floor(max(log.Quantiles)+1),ceiling((max(log.Quantiles)+1)/10)), cex.axis=cex.axis)
+				axis(2, at=seq(0,floor(YlimMax+1),ceiling((YlimMax+1)/10)), labels=seq(0,floor((YlimMax+1)),ceiling((YlimMax+1)/10)), cex.axis=cex.axis)
+				
+				# plot the confidence interval of QQ-plot
+				if(conf.int)	polygon(c(log.Quantiles[index],log.Quantiles),c(-log10(c05)[index],-log10(c95)),col=conf.int.col,border=conf.int.col)
+				
+				for(i in 1:R){
+					if(verbose)	print(paste("Multraits_QQ Plotting ",taxa[i],"...",sep=""))
+					P.values=as.numeric(Pmap[,i+2])
+					P.values=P.values[!is.na(P.values)]
+					if(LOG10){
+						P.values=P.values[P.values>=0]
+						P.values=P.values[P.values<=1]
+						N=length(P.values)
+						P.values=P.values[order(P.values)]
+					}else{
+						N=length(P.values)
+						P.values=P.values[order(P.values,decreasing=TRUE)]
+					}
+					if(LOG10){
+						log.P.values <- -log10(P.values)
+					}else{
+						log.P.values <- P.values
+					}
+
+						
+					if((i == 1) & !is.null(threshold.col)){par(xpd=FALSE); abline(a = 0, b = 1, col = threshold.col[1],lwd=2); par(xpd=TRUE)}
+					points(log.Quantiles, log.P.values, col = t(col)[i],pch=19,cex=cex[3])
+						
+					if(!is.null(threshold)){
+						if(sum(threshold!=0)==length(threshold)){
+							thre.line=-log10(min(threshold))
+							if(amplify==TRUE){
+								thre.index=which(log.P.values>=thre.line)
+								if(length(thre.index)!=0){
+								
+									# cover the points that exceed the threshold with the color "white"
+									points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,cex=cex[3])
+									if(is.null(signal.col)){
+										points(log.Quantiles[thre.index],log.P.values[thre.index],col = t(col)[i],pch=signal.pch[1],cex=signal.cex[1])
+									}else{
+										points(log.Quantiles[thre.index],log.P.values[thre.index],col = signal.col[1],pch=signal.pch[1],cex=signal.cex[1])
+									}
+								}
+							}
+						}
+					}
+				}
+				if(box)	box()
+				if(file.output) dev.off()
+			}
+		}else{
+			for(i in 1:R){
+				if(verbose)	print(paste("Q_Q Plotting ",taxa[i],"...",sep=""))
+				if(file.output){
+					if(file=="jpg")	jpeg(paste("QQplot.",taxa[i],".jpg",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
+					if(file=="pdf")	pdf(paste("QQplot.",taxa[i],".pdf",sep=""), width = 5.5,height=5.5)
+					if(file=="tiff")	tiff(paste("QQplot.",taxa[i],".tiff",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi)
+					par(mar = c(5,5,4,2),xpd=TRUE)
+				}else{
+					if(is.null(dev.list()))	dev.new(width = 5.5, height = 5.5)
+					par(xpd=TRUE)
+				}
+				P.values=as.numeric(Pmap[,i+2])
+				P.values=P.values[!is.na(P.values)]
+				if(LOG10){
+					P.values=P.values[P.values>0]
+					P.values=P.values[P.values<=1]
+					N=length(P.values)
+					P.values=P.values[order(P.values)]
+				}else{
+					N=length(P.values)
+					P.values=P.values[order(P.values,decreasing=TRUE)]
+				}
+				p_value_quantiles=(1:length(P.values))/(length(P.values)+1)
+				log.Quantiles <- -log10(p_value_quantiles)
+				if(LOG10){
+					log.P.values <- -log10(P.values)
+				}else{
+					log.P.values <- P.values
+				}
+				
+				#calculate the confidence interval of QQ-plot
+				if(conf.int){
+					N1=length(log.Quantiles)
+					c95 <- rep(NA,N1)
+					c05 <- rep(NA,N1)
+					for(j in 1:N1){
+						xi=ceiling((10^-log.Quantiles[j])*N)
+						if(xi==0)xi=1
+						c95[j] <- qbeta(0.95,xi,N-xi+1)
+						c05[j] <- qbeta(0.05,xi,N-xi+1)
+					}
+					index=length(c95):1
+				}else{
+					c05 <- 1
+					c95 <- 1
+				}
+				YlimMax <- max(floor(max(max(-log10(c05)), max(-log10(c95)))+1), floor(max(log.P.values)+1))
+				plot(NULL, xlim = c(0,floor(max(log.Quantiles)+1)), axes=FALSE, cex.axis=cex.axis, cex.lab=1.2,ylim=c(0,YlimMax),xlab =expression(Expected~~-log[10](italic(p))), ylab = expression(Observed~~-log[10](italic(p))), main = paste("QQplot of",taxa[i]))
+				axis(1, at=seq(0,floor(max(log.Quantiles)+1),ceiling((max(log.Quantiles)+1)/10)), labels=seq(0,floor(max(log.Quantiles)+1),ceiling((max(log.Quantiles)+1)/10)), cex.axis=cex.axis)
+				axis(2, at=seq(0,YlimMax,ceiling(YlimMax/10)), labels=seq(0,YlimMax,ceiling(YlimMax/10)), cex.axis=cex.axis)
+				
+				#plot the confidence interval of QQ-plot
+				if(conf.int)	polygon(c(log.Quantiles[index],log.Quantiles),c(-log10(c05)[index],-log10(c95)),col=conf.int.col,border=conf.int.col)
+				
+				if(!is.null(threshold.col)){par(xpd=FALSE); abline(a = 0, b = 1, col = threshold.col[1],lwd=2); par(xpd=TRUE)}
+				points(log.Quantiles, log.P.values, col = col[1],pch=19,cex=cex[3])
+				
+				if(!is.null(threshold)){
+					if(sum(threshold!=0)==length(threshold)){
+						thre.line=-log10(min(threshold))
+						if(amplify==TRUE){
+							thre.index=which(log.P.values>=thre.line)
+							if(length(thre.index)!=0){
+							
+								#cover the points that exceed the threshold with the color "white"
+								points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,cex=cex[3])
+								if(is.null(signal.col)){
+									points(log.Quantiles[thre.index],log.P.values[thre.index],col = col[1],pch=signal.pch[1],cex=signal.cex[1])
+								}else{
+									points(log.Quantiles[thre.index],log.P.values[thre.index],col = signal.col[1],pch=signal.pch[1],cex=signal.cex[1])
+								}
+							}
+						}
+					}
+				}
+				if(box)	box()
+				if(file.output) dev.off()
+			}
+		}
+	}
+	if(file.output & verbose)	print(paste("Plots are stored in: ", getwd(), sep=""))
+}