#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)