Datasets
Models
Diff of
/modas/utils/CMplot.r
[000000]
..
[a43cea]
Switch to side-by-side view
--- 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="")) +}