R/GAPIT.Circle.Manhattan.Plot.R
In jiabowang/GAPIT3: GAPIT (Genomic Association and Prediction Integrated Tool)
Defines functions
GAPIT.Circle.Manhattan.Plot
Densitplot
circle.plot
########## These three functions come from MVP package, Jiabo did some modifications
########## Following Apache License, we thank MVP developper to build these functions.
########## 1 creat P value scale in addtitional chromsome
########## 2 set col same as GAPIT
########## 3
circle.plot <- function(myr,type="l",x=NULL,lty=1,lwd=1,col="black",add=TRUE,n.point=1000)
{
graphics::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)
graphics::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("darkblue", "white", "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
)
{ #print(head(map))
map <- as.matrix(map)
map <- map[!is.na(map[, 2]), ]
map <- map[!is.na(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)
#print(chr.num)
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
#print(chr.num)
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){
col.index[[i]] = 1
}else{
cut.r <- cut(pos.x[[i]], cut.len, labels=FALSE)
eachbin.num <- table(cut.r)
#print(eachbin.num)
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
}
#print(col)
#print(maxbin.num)
col = grDevices::colorRampPalette(col)(maxbin.num)
col.seg=NULL
for(i in 1 : length(chr.num)){
if(plot) graphics::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) graphics::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) graphics::mtext(at=seq(band, length(chr.num) * band, band),text=paste("Chr", chr.num, sep=""), side=2, las=2, font=1, cex=0.6, line=0.2)
if(plot) graphics::axis(3, at=seq(0, chorm.maxlen, length=10), labels=c(NA, paste(round((seq(0, chorm.maxlen, length=10))[-1] / 1e6, 0), "Mb", sep="")),
font=1, 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)
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) graphics::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))
}
GAPIT.Circle.Manhattan.Plot <- function(
Pmap,
col=c("#377EB8", "#4DAF4A", "#984EA3", "#FF7F00"),
#col=c("darkgreen", "darkblue", "darkyellow", "darkred"),
bin.size=1e6,
bin.max=NULL,
pch=19,
band=1,
cir.band=0.5,
H=1.5,
ylim=NULL,
cex.axis=1,
plot.type="c",
multracks=TRUE,
cex=c(0.5,0.8,1),
r=0.3,
xlab="Chromosome",
ylab=expression(-log[10](italic(p))),
xaxs="i",
yaxs="r",
outward=TRUE,
threshold = 0.01,
threshold.col="red",
threshold.lwd=1,
threshold.lty=2,
amplify= TRUE, # is that available for remark signal pch col
chr.labels=NULL,
signal.cex = 2,
signal.pch = 8,
signal.col="red",
signal.line=NULL,
cir.chr=TRUE,
cir.chr.h=1.3,
chr.den.col=c("darkgreen", "yellow", "red"),
#chr.den.col=c(126,177,153),
cir.legend=TRUE,
cir.legend.cex=0.8,
cir.legend.col="grey45",
LOG10=TRUE,
box=FALSE,
conf.int.col="grey",
file.output=TRUE,
file="pdf",
dpi=300,
xz=NULL,
memo=""
)
{ #print("Starting Circular-Manhattan plot!",quote=F)
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="")
col=rep(c( '#FF6A6A', '#FAC863', '#99C794', '#6699CC', '#C594C5'),ceiling(length(taxa)/5))
legend.bit=round(nrow(Pmap)/30)
numeric.chr <- as.numeric(Pmap[, 1])
options(warn = 0)
max.chr <- max(numeric.chr, na.rm=TRUE)
aa=Pmap[1:legend.bit,]
aa[,2]=max.chr+1
#print(aa[,3])
aa[,3]=sample(1:10^7.5,legend.bit)
aa[,-c(1:3)]=0
Pmap=rbind(Pmap,aa)
#print(unique(Pmap[,2]))
#SNP-Density plot
if("d" %in% plot.type){
print("SNP_Density Plotting...")
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("SNP_Density.",paste(taxa,collapse="."),".jpg",sep=""), width = 9*dpi,height=7*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.SNP_Density", taxa,".pdf" ,sep=""), width = 9,height=7)
if(file=="tiff") grDevices::tiff(paste("SNP_Density.",paste(taxa,collapse="."),".tiff",sep=""), width = 9*dpi,height=7*dpi,res=dpi)
graphics::par(xpd=TRUE)
}else{
if(is.null(grDevices::dev.list())) grDevices::dev.new(width = 9,height=7)
graphics::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) grDevices::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]
Pmap[is.na(Pmap)]=1
#print(dim(Pmap))
#scale and adjust the parameters
cir.chr.h <- cir.chr.h/5
cir.band <- cir.band/5
threshold=threshold/nrow(Pmap)
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))
#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)])
#print(dim(pvalueT))
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)))
#print(col)
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
#print(Nchr)
#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{
pvalueT <- abs(pvalueT)
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
}
#print(dim(pvalueT))
if(is.null(xz)){
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.lty=rep(2,length(signal.line.index))
}else{
signal.line.index=as.numeric(as.vector(xz[,1]))
signal.lty=as.numeric(as.vector(xz[,2]))
}#end is.null(xz)
signal.line.index <- pvalue.posN[signal.line.index]
}
if("c" %in% plot.type)
{
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("GAPIT.Manhattan.Multiple.Plot.circular.jpg",sep=""), width = 8*dpi,height=8*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.Manhattans_Circular.pdf" ,sep=""), width = 10,height=10)
if(file=="tiff") grDevices::tiff(paste("GAPIT.Manhattan.Multiple.Plot.circular.tiff",sep=""), width = 8*dpi,height=8*dpi,res=dpi)
}
if(!file.output){
if(!is.null(grDevices::dev.list())) grDevices::dev.new(width=8, height=8)
graphics::par(pty="s", xpd=TRUE, mar=c(1,1,1,1))
}
graphics::par(pty="s", xpd=TRUE, mar=c(1,1,1,1))
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)
#print(signal.line)
#print(dim(pvalueT))
#print(head(pvalueT))
#print(dim(xz))
#print(xz)
#print(head(pvalue.posN))
graphics::segments(X1chr,Y1chr,X2chr,Y2chr,lty=signal.lty,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)
#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])))
}else{
Max <- ceiling(max(pvalue[pvalue!=Inf]))
if(Max<=1)
Max <- max(pvalue[pvalue!=Inf])
}
}else{
Max <- ylim[2]
}
Cpvalue <- (H*logpvalue/Max)
if(outward==TRUE){
if(cir.chr==TRUE){
#plot the boundary which represents the chromosomes
polygon.num <- 1000
#print(length(chr))
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)
#print(length(X1chr))
if(is.null(chr.den.col)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::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)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
}
}
}
}
if(cir.density){
graphics::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
)
graphics::legend(
x=RR+4*cir.chr.h,
y=(RR+4*cir.chr.h)/2,
horiz=F,
title="Density", 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)
}
}
#plot the y axis of legend for each trait
if(cir.legend==TRUE){
#try to get the number after radix point
if(Max<=1) {
round.n=nchar(as.character(10^(-ceiling(-log10(Max)))))-1
}else{
round.n=1
}
graphics::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)
graphics::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')
graphics::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')
graphics::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')
graphics::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')
graphics::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.75)+cir.band*(i-1),round(Max*0.25,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
graphics::text(-r/15,r+H*(i-0.5)+cir.band*(i-1),round(Max*0.5,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
graphics::text(-r/15,r+H*(i-0.25)+cir.band*(i-1),round(Max*0.75,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(-r/15,r+H*(i-0)+cir.band*(i-1),round(Max*1,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(r/5,0.4*(i-1),taxa[i],adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
}
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)
# plot point in figure
graphics::points(X[1:(length(X)-legend.bit)],Y[1:(length(Y)-legend.bit)],pch=19,cex=cex[1],col=rep(rep(colx,N[i]),add[[i]]))
# plot significant line
if(!is.null(threshold)){
if(sum(threshold!=0)==length(threshold)){
for(thr in 1:length(threshold)){
significantline1=ifelse(LOG10, H*(-log10(threshold[thr]))/Max, H*(threshold[thr])/Max)
#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)
# plot significant line
if(significantline1<H){
#lines(s1X,s1Y,type="l",col=threshold.col,lwd=threshold.col,lty=threshold.lty)
#if(thr==length(threshold))circle.plot(myr=(significantline1+r+H*(i-1)+cir.band*(i-1)),col="black",lwd=threshold.lwd[thr],lty=threshold.lty[thr])
#print("!!!!!")
circle.plot(myr=(significantline1+r+H*(i-1)+cir.band*(i-1)),col=threshold.col[thr],lwd=threshold.lwd[thr],lty=threshold.lty[thr])
#circle.plot(myr=(significantline1+r+H*(i-1)+cir.band*(i-1)),col="black",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)))/Max
}else{
threshold <- sort(threshold, decreasing=TRUE)
significantline1=H*(min(threshold))/Max
}
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"
graphics::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]))/Max
}else{
significantline1=H*(threshold[ll])/Max
}
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]))/Max
significantline1=H*(-log10(threshold[ll]))/Max
}else{
significantline0=H*(threshold[ll-1])/Max
significantline1=H*(threshold[ll])/Max
}
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)){
# print(signal.pch)
graphics::points(HX1,HY1,pch=signal.pch,cex=signal.cex[ll]*cex[1],col=rep(rep(colx,N[i]),add[[i]])[p_amp.index])
}else{
# print(signal.pch)
graphics::points(HX1,HY1,pch=signal.pch,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)){
#print(length(ticks))
for(i in 1:(length(ticks)-1)){
angle=360*(1-(ticks-round(band/2))[i]/TotalN)
graphics::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)
graphics::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)
graphics::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)
graphics::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)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::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)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
}
}
}
}
if(cir.density){
graphics::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
)
graphics::legend(
x=RR+4*cir.chr.h,
y=(RR+4*cir.chr.h)/2,
title="Density", 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)
}
}
if(cir.legend==TRUE){
#try to get the number after radix point
if(Max<=1) {
round.n=nchar(as.character(10^(-ceiling(-log10(Max)))))-1
}else{
round.n=2
}
graphics::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)
graphics::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')
graphics::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')
graphics::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')
graphics::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')
graphics::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')
graphics::text(-r/15,r+H*(i-0.25)+cir.band*(i-1),round(Max*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(Max*0.5,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
graphics::text(-r/15,r+H*(i-0.75)+cir.band*(i-1),round(Max*0.75,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(-r/15,r+H*(i-1)+cir.band*(i-1),round(Max*1,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(r,0.4*(i-1),taxa[i],adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
}
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]]))
graphics::points(X[1:(length(X)-legend.bit)],Y[1:(length(Y)-legend.bit)],pch=19,cex=cex[1],col=rep(rep(colx,N[i]),add[[i]]))
if(!is.null(threshold)){
if(sum(threshold!=0)==length(threshold)){
for(thr in 1:length(threshold)){
significantline1=ifelse(LOG10, H*(-log10(threshold[thr]))/Max, H*(threshold[thr])/Max)
#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)))/Max
}else{
threshold <- sort(threshold, decreasing=TRUE)
significantline1=H*(min(threshold))/Max
}
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"
graphics::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]))/Max
}else{
significantline1=H*(threshold[ll])/Max
}
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]))/Max
significantline1=H*(-log10(threshold[ll]))/Max
}else{
significantline0=H*(threshold[ll-1])/Max
significantline1=H*(threshold[ll])/Max
}
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)){
graphics::points(HX1,HY1,pch=signal.pch,cex=signal.cex[ll]*cex[1],col=rep(rep(colx,N[i]),add[[i]])[p_amp.index])
}else{
graphics::points(HX1,HY1,pch=signal.pch,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)-1)){
angle=360*(1-(ticks-round(band/2))[i]/TotalN)
graphics::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)
graphics::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)
graphics::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)
graphics::text(ticks1[i],ticks2[i],chr.labels[i],srt=angle,font=2,cex=cex.axis)
}
}
}
}
}
taxa=append("Centre",taxa,)
taxa_col=rep("black",R)
taxa_col=append("red",taxa_col)
for(j in 1:(R+1)){
graphics::text(r/5,0.4*(j-1),taxa[j],adj=1,col=taxa_col[j],cex=cir.legend.cex,font=2)
}
taxa=taxa[-1]
if(file.output) grDevices::dev.off()
}
if("q" %in% plot.type){
#print("Starting QQ-plot!",quote=F)
amplify=FALSE
if(multracks){
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("GAPIT.Multracks.QQ.plot.jpg",sep=""), width = R*2.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.QQs_Tracks.pdf",sep=""), width = R*2.5,height=5.5)
if(file=="tiff") grDevices::tiff(paste("GAPIT.Multracks.QQ.plot.tiff",sep=""), width = R*2.5*dpi,height=5.5*dpi,res=dpi)
graphics::par(mfcol=c(1,R),mar = c(0,1,4,1.5),oma=c(3,5,0,0),xpd=TRUE)
}else{
if(is.null(grDevices::dev.list())) grDevices::dev.new(width = 2.5*R, height = 5.5)
graphics::par(xpd=TRUE)
}
for(i in 1:R){
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] <- stats::qbeta(0.95,xi,N-xi+1)
c05[j] <- stats::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])
graphics::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)
graphics::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) graphics::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)){
graphics::par(xpd=FALSE);
graphics::abline(a = 0, b = 1, col = threshold.col[1],lwd=2);
graphics::par(xpd=TRUE)
}
graphics::points(log.Quantiles, log.P.values, col = col[1],pch=1,cex=cex[3])
#print(max(log.Quantiles))
# print(length(log.Quantiles))
# print(length(log.P.values))
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"
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,cex=cex[3])
if(is.null(signal.col)){
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = col[1],pch=signal.pch[1],cex=signal.cex[1])
}else{
graphics::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) grDevices::dev.off()
if(R > 1){
#qq_col=rainbow(R)
qq_col=rep(c( '#FF6A6A', '#FAC863', '#99C794', '#6699CC', '#C594C5'),ceiling(R/5))
signal.col <- NULL
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("GAPIT.Multiple.QQ.plot.symphysic.jpg",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.QQs_Symphysic.pdf",sep=""), width = 5.5,height=5.5)
if(file=="tiff") grDevices::tiff(paste("GAPIT.Multiple.QQ.plot.symphysic.tiff",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi)
graphics::par(mar = c(5,5,4,2),xpd=TRUE)
}else{
grDevices::dev.new(width = 5.5, height = 5.5)
graphics::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)
# 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] <- stats::qbeta(0.95,xi,N-xi+1)
c05[j] <- stats::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 = "QQ plot")
#legend("topleft",taxa,col=t(col)[1:R],pch=1,pt.lwd=2,text.font=6,box.col=NA)
graphics::legend("topleft",taxa,col=qq_col[1:R],pch=1,pt.lwd=3,text.font=6,box.col=NA)
graphics::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)
graphics::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)
#print(log.Quantiles[index])
#print(index)
#print(length(log.Quantiles))
# plot the confidence interval of QQ-plot
if(conf.int) graphics::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){
#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)]
}
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
}
if((i == 1) & !is.null(threshold.col)){
graphics::par(xpd=FALSE);
graphics::abline(a = 0, b = 1, col = threshold.col[1],lwd=2);
graphics::par(xpd=TRUE)}
#print(length(log.Quantiles))
#print("!!!!!")
#points(log.Quantiles, log.P.values, col = t(col)[i],pch=1,lwd=3,cex=cex[3])
graphics::points(log.Quantiles, log.P.values, col = qq_col[i],pch=1,lwd=3,cex=cex[3])
#print(max(log.Quantiles))
#
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"
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,lwd=3,cex=cex[3])
if(is.null(signal.col)){
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = t(col)[i],pch=signal.pch[1],cex=signal.cex[1])
}else{
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = signal.col[1],pch=signal.pch[1],cex=signal.cex[1])
}
}
}
}
}
}
box()
if(file.output) grDevices::dev.off()
}
}else{
for(i in 1:R){
print(paste("Q_Q Plotting ",taxa[i],"...",sep=""))
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("QQplot.",taxa[i],".jpg",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.Q_Q.",taxa[i],".pdf",sep=""), width = 5.5,height=5.5)
if(file=="tiff") grDevices::tiff(paste("QQplot.",taxa[i],".tiff",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi)
graphics::par(mar = c(5,5,4,2),xpd=TRUE)
}else{
if(is.null(grDevices::dev.list())) grDevices::dev.new(width = 5.5, height = 5.5)
graphics::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] <- stats::qbeta(0.95,xi,N-xi+1)
c05[j] <- stats::qbeta(0.05,xi,N-xi+1)
}
index=length(c95):1
}else{
c05 <- 1
c95 <- 1
}
#print(max(log.Quantiles))
#print("@@@@@")
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]))
graphics::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)
graphics::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
#print(log.Quantiles[index])
qq_col = grDevices::rainbow(R)
#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(conf.int) graphics::polygon(c(log.Quantiles[index],log.Quantiles),c(-log10(c05)[index],-log10(c95)),col=qq_col[i],border=conf.int.col)
if(!is.null(threshold.col)){
graphics::par(xpd=FALSE);
graphics::abline(a = 0, b = 1, col = threshold.col[1],lwd=2);
graphics::par(xpd=TRUE)
}
graphics::points(log.Quantiles, log.P.values, col = col[1],pch=19,cex=2)
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){
#print("!!!!")
#cover the points that exceed the threshold with the color "white"
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,lwd=3,cex=cex[3])
if(is.null(signal.col)){
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = col[1],pch=signal.pch[1],cex=signal.cex[1])
}else{
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = signal.col[1],pch=signal.pch[1],cex=signal.cex[1])
}
}
}
}
}
box()
if(file.output) grDevices::dev.off()
}
}
print("Multiple QQ plot has been finished!")
}
}#End of Whole function
#}
jiabowang/GAPIT3 documentation built on March 6, 2025, 2:21 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions GAPIT.Circle.Manhattan.Plot Densitplot circle.plot
########## These three functions come from MVP package, Jiabo did some modifications
########## Following Apache License, we thank MVP developper to build these functions.
########## 1 creat P value scale in addtitional chromsome
########## 2 set col same as GAPIT
########## 3
circle.plot <- function(myr,type="l",x=NULL,lty=1,lwd=1,col="black",add=TRUE,n.point=1000)
{
graphics::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)
graphics::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("darkblue", "white", "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
)
{ #print(head(map))
map <- as.matrix(map)
map <- map[!is.na(map[, 2]), ]
map <- map[!is.na(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)
#print(chr.num)
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
#print(chr.num)
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){
col.index[[i]] = 1
}else{
cut.r <- cut(pos.x[[i]], cut.len, labels=FALSE)
eachbin.num <- table(cut.r)
#print(eachbin.num)
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
}
#print(col)
#print(maxbin.num)
col = grDevices::colorRampPalette(col)(maxbin.num)
col.seg=NULL
for(i in 1 : length(chr.num)){
if(plot) graphics::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) graphics::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) graphics::mtext(at=seq(band, length(chr.num) * band, band),text=paste("Chr", chr.num, sep=""), side=2, las=2, font=1, cex=0.6, line=0.2)
if(plot) graphics::axis(3, at=seq(0, chorm.maxlen, length=10), labels=c(NA, paste(round((seq(0, chorm.maxlen, length=10))[-1] / 1e6, 0), "Mb", sep="")),
font=1, 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)
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) graphics::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))
}
GAPIT.Circle.Manhattan.Plot <- function(
Pmap,
col=c("#377EB8", "#4DAF4A", "#984EA3", "#FF7F00"),
#col=c("darkgreen", "darkblue", "darkyellow", "darkred"),
bin.size=1e6,
bin.max=NULL,
pch=19,
band=1,
cir.band=0.5,
H=1.5,
ylim=NULL,
cex.axis=1,
plot.type="c",
multracks=TRUE,
cex=c(0.5,0.8,1),
r=0.3,
xlab="Chromosome",
ylab=expression(-log[10](italic(p))),
xaxs="i",
yaxs="r",
outward=TRUE,
threshold = 0.01,
threshold.col="red",
threshold.lwd=1,
threshold.lty=2,
amplify= TRUE, # is that available for remark signal pch col
chr.labels=NULL,
signal.cex = 2,
signal.pch = 8,
signal.col="red",
signal.line=NULL,
cir.chr=TRUE,
cir.chr.h=1.3,
chr.den.col=c("darkgreen", "yellow", "red"),
#chr.den.col=c(126,177,153),
cir.legend=TRUE,
cir.legend.cex=0.8,
cir.legend.col="grey45",
LOG10=TRUE,
box=FALSE,
conf.int.col="grey",
file.output=TRUE,
file="pdf",
dpi=300,
xz=NULL,
memo=""
)
{ #print("Starting Circular-Manhattan plot!",quote=F)
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="")
col=rep(c( '#FF6A6A', '#FAC863', '#99C794', '#6699CC', '#C594C5'),ceiling(length(taxa)/5))
legend.bit=round(nrow(Pmap)/30)
numeric.chr <- as.numeric(Pmap[, 1])
options(warn = 0)
max.chr <- max(numeric.chr, na.rm=TRUE)
aa=Pmap[1:legend.bit,]
aa[,2]=max.chr+1
#print(aa[,3])
aa[,3]=sample(1:10^7.5,legend.bit)
aa[,-c(1:3)]=0
Pmap=rbind(Pmap,aa)
#print(unique(Pmap[,2]))
#SNP-Density plot
if("d" %in% plot.type){
print("SNP_Density Plotting...")
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("SNP_Density.",paste(taxa,collapse="."),".jpg",sep=""), width = 9*dpi,height=7*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.SNP_Density", taxa,".pdf" ,sep=""), width = 9,height=7)
if(file=="tiff") grDevices::tiff(paste("SNP_Density.",paste(taxa,collapse="."),".tiff",sep=""), width = 9*dpi,height=7*dpi,res=dpi)
graphics::par(xpd=TRUE)
}else{
if(is.null(grDevices::dev.list())) grDevices::dev.new(width = 9,height=7)
graphics::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) grDevices::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]
Pmap[is.na(Pmap)]=1
#print(dim(Pmap))
#scale and adjust the parameters
cir.chr.h <- cir.chr.h/5
cir.band <- cir.band/5
threshold=threshold/nrow(Pmap)
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))
#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)])
#print(dim(pvalueT))
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)))
#print(col)
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
#print(Nchr)
#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{
pvalueT <- abs(pvalueT)
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
}
#print(dim(pvalueT))
if(is.null(xz)){
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.lty=rep(2,length(signal.line.index))
}else{
signal.line.index=as.numeric(as.vector(xz[,1]))
signal.lty=as.numeric(as.vector(xz[,2]))
}#end is.null(xz)
signal.line.index <- pvalue.posN[signal.line.index]
}
if("c" %in% plot.type)
{
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("GAPIT.Manhattan.Multiple.Plot.circular.jpg",sep=""), width = 8*dpi,height=8*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.Manhattans_Circular.pdf" ,sep=""), width = 10,height=10)
if(file=="tiff") grDevices::tiff(paste("GAPIT.Manhattan.Multiple.Plot.circular.tiff",sep=""), width = 8*dpi,height=8*dpi,res=dpi)
}
if(!file.output){
if(!is.null(grDevices::dev.list())) grDevices::dev.new(width=8, height=8)
graphics::par(pty="s", xpd=TRUE, mar=c(1,1,1,1))
}
graphics::par(pty="s", xpd=TRUE, mar=c(1,1,1,1))
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)
#print(signal.line)
#print(dim(pvalueT))
#print(head(pvalueT))
#print(dim(xz))
#print(xz)
#print(head(pvalue.posN))
graphics::segments(X1chr,Y1chr,X2chr,Y2chr,lty=signal.lty,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)
#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])))
}else{
Max <- ceiling(max(pvalue[pvalue!=Inf]))
if(Max<=1)
Max <- max(pvalue[pvalue!=Inf])
}
}else{
Max <- ylim[2]
}
Cpvalue <- (H*logpvalue/Max)
if(outward==TRUE){
if(cir.chr==TRUE){
#plot the boundary which represents the chromosomes
polygon.num <- 1000
#print(length(chr))
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)
#print(length(X1chr))
if(is.null(chr.den.col)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::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)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
}
}
}
}
if(cir.density){
graphics::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
)
graphics::legend(
x=RR+4*cir.chr.h,
y=(RR+4*cir.chr.h)/2,
horiz=F,
title="Density", 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)
}
}
#plot the y axis of legend for each trait
if(cir.legend==TRUE){
#try to get the number after radix point
if(Max<=1) {
round.n=nchar(as.character(10^(-ceiling(-log10(Max)))))-1
}else{
round.n=1
}
graphics::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)
graphics::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')
graphics::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')
graphics::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')
graphics::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')
graphics::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.75)+cir.band*(i-1),round(Max*0.25,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
graphics::text(-r/15,r+H*(i-0.5)+cir.band*(i-1),round(Max*0.5,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
graphics::text(-r/15,r+H*(i-0.25)+cir.band*(i-1),round(Max*0.75,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(-r/15,r+H*(i-0)+cir.band*(i-1),round(Max*1,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(r/5,0.4*(i-1),taxa[i],adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
}
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)
# plot point in figure
graphics::points(X[1:(length(X)-legend.bit)],Y[1:(length(Y)-legend.bit)],pch=19,cex=cex[1],col=rep(rep(colx,N[i]),add[[i]]))
# plot significant line
if(!is.null(threshold)){
if(sum(threshold!=0)==length(threshold)){
for(thr in 1:length(threshold)){
significantline1=ifelse(LOG10, H*(-log10(threshold[thr]))/Max, H*(threshold[thr])/Max)
#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)
# plot significant line
if(significantline1<H){
#lines(s1X,s1Y,type="l",col=threshold.col,lwd=threshold.col,lty=threshold.lty)
#if(thr==length(threshold))circle.plot(myr=(significantline1+r+H*(i-1)+cir.band*(i-1)),col="black",lwd=threshold.lwd[thr],lty=threshold.lty[thr])
#print("!!!!!")
circle.plot(myr=(significantline1+r+H*(i-1)+cir.band*(i-1)),col=threshold.col[thr],lwd=threshold.lwd[thr],lty=threshold.lty[thr])
#circle.plot(myr=(significantline1+r+H*(i-1)+cir.band*(i-1)),col="black",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)))/Max
}else{
threshold <- sort(threshold, decreasing=TRUE)
significantline1=H*(min(threshold))/Max
}
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"
graphics::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]))/Max
}else{
significantline1=H*(threshold[ll])/Max
}
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]))/Max
significantline1=H*(-log10(threshold[ll]))/Max
}else{
significantline0=H*(threshold[ll-1])/Max
significantline1=H*(threshold[ll])/Max
}
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)){
# print(signal.pch)
graphics::points(HX1,HY1,pch=signal.pch,cex=signal.cex[ll]*cex[1],col=rep(rep(colx,N[i]),add[[i]])[p_amp.index])
}else{
# print(signal.pch)
graphics::points(HX1,HY1,pch=signal.pch,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)){
#print(length(ticks))
for(i in 1:(length(ticks)-1)){
angle=360*(1-(ticks-round(band/2))[i]/TotalN)
graphics::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)
graphics::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)
graphics::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)
graphics::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)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::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)){
graphics::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){
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col="grey",border="grey")
}else{
graphics::polygon(c(rev(X1chr),X2chr),c(rev(Y1chr),Y2chr),col=chr.den.col,border=chr.den.col)
}
}
}
}
if(cir.density){
graphics::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
)
graphics::legend(
x=RR+4*cir.chr.h,
y=(RR+4*cir.chr.h)/2,
title="Density", 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)
}
}
if(cir.legend==TRUE){
#try to get the number after radix point
if(Max<=1) {
round.n=nchar(as.character(10^(-ceiling(-log10(Max)))))-1
}else{
round.n=2
}
graphics::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)
graphics::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')
graphics::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')
graphics::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')
graphics::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')
graphics::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')
graphics::text(-r/15,r+H*(i-0.25)+cir.band*(i-1),round(Max*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(Max*0.5,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
graphics::text(-r/15,r+H*(i-0.75)+cir.band*(i-1),round(Max*0.75,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(-r/15,r+H*(i-1)+cir.band*(i-1),round(Max*1,round.n),adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
#text(r,0.4*(i-1),taxa[i],adj=1,col=cir.legend.col,cex=cir.legend.cex,font=2)
}
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]]))
graphics::points(X[1:(length(X)-legend.bit)],Y[1:(length(Y)-legend.bit)],pch=19,cex=cex[1],col=rep(rep(colx,N[i]),add[[i]]))
if(!is.null(threshold)){
if(sum(threshold!=0)==length(threshold)){
for(thr in 1:length(threshold)){
significantline1=ifelse(LOG10, H*(-log10(threshold[thr]))/Max, H*(threshold[thr])/Max)
#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)))/Max
}else{
threshold <- sort(threshold, decreasing=TRUE)
significantline1=H*(min(threshold))/Max
}
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"
graphics::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]))/Max
}else{
significantline1=H*(threshold[ll])/Max
}
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]))/Max
significantline1=H*(-log10(threshold[ll]))/Max
}else{
significantline0=H*(threshold[ll-1])/Max
significantline1=H*(threshold[ll])/Max
}
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)){
graphics::points(HX1,HY1,pch=signal.pch,cex=signal.cex[ll]*cex[1],col=rep(rep(colx,N[i]),add[[i]])[p_amp.index])
}else{
graphics::points(HX1,HY1,pch=signal.pch,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)-1)){
angle=360*(1-(ticks-round(band/2))[i]/TotalN)
graphics::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)
graphics::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)
graphics::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)
graphics::text(ticks1[i],ticks2[i],chr.labels[i],srt=angle,font=2,cex=cex.axis)
}
}
}
}
}
taxa=append("Centre",taxa,)
taxa_col=rep("black",R)
taxa_col=append("red",taxa_col)
for(j in 1:(R+1)){
graphics::text(r/5,0.4*(j-1),taxa[j],adj=1,col=taxa_col[j],cex=cir.legend.cex,font=2)
}
taxa=taxa[-1]
if(file.output) grDevices::dev.off()
}
if("q" %in% plot.type){
#print("Starting QQ-plot!",quote=F)
amplify=FALSE
if(multracks){
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("GAPIT.Multracks.QQ.plot.jpg",sep=""), width = R*2.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.QQs_Tracks.pdf",sep=""), width = R*2.5,height=5.5)
if(file=="tiff") grDevices::tiff(paste("GAPIT.Multracks.QQ.plot.tiff",sep=""), width = R*2.5*dpi,height=5.5*dpi,res=dpi)
graphics::par(mfcol=c(1,R),mar = c(0,1,4,1.5),oma=c(3,5,0,0),xpd=TRUE)
}else{
if(is.null(grDevices::dev.list())) grDevices::dev.new(width = 2.5*R, height = 5.5)
graphics::par(xpd=TRUE)
}
for(i in 1:R){
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] <- stats::qbeta(0.95,xi,N-xi+1)
c05[j] <- stats::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])
graphics::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)
graphics::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) graphics::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)){
graphics::par(xpd=FALSE);
graphics::abline(a = 0, b = 1, col = threshold.col[1],lwd=2);
graphics::par(xpd=TRUE)
}
graphics::points(log.Quantiles, log.P.values, col = col[1],pch=1,cex=cex[3])
#print(max(log.Quantiles))
# print(length(log.Quantiles))
# print(length(log.P.values))
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"
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,cex=cex[3])
if(is.null(signal.col)){
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = col[1],pch=signal.pch[1],cex=signal.cex[1])
}else{
graphics::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) grDevices::dev.off()
if(R > 1){
#qq_col=rainbow(R)
qq_col=rep(c( '#FF6A6A', '#FAC863', '#99C794', '#6699CC', '#C594C5'),ceiling(R/5))
signal.col <- NULL
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("GAPIT.Multiple.QQ.plot.symphysic.jpg",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.QQs_Symphysic.pdf",sep=""), width = 5.5,height=5.5)
if(file=="tiff") grDevices::tiff(paste("GAPIT.Multiple.QQ.plot.symphysic.tiff",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi)
graphics::par(mar = c(5,5,4,2),xpd=TRUE)
}else{
grDevices::dev.new(width = 5.5, height = 5.5)
graphics::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)
# 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] <- stats::qbeta(0.95,xi,N-xi+1)
c05[j] <- stats::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 = "QQ plot")
#legend("topleft",taxa,col=t(col)[1:R],pch=1,pt.lwd=2,text.font=6,box.col=NA)
graphics::legend("topleft",taxa,col=qq_col[1:R],pch=1,pt.lwd=3,text.font=6,box.col=NA)
graphics::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)
graphics::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)
#print(log.Quantiles[index])
#print(index)
#print(length(log.Quantiles))
# plot the confidence interval of QQ-plot
if(conf.int) graphics::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){
#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)]
}
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
}
if((i == 1) & !is.null(threshold.col)){
graphics::par(xpd=FALSE);
graphics::abline(a = 0, b = 1, col = threshold.col[1],lwd=2);
graphics::par(xpd=TRUE)}
#print(length(log.Quantiles))
#print("!!!!!")
#points(log.Quantiles, log.P.values, col = t(col)[i],pch=1,lwd=3,cex=cex[3])
graphics::points(log.Quantiles, log.P.values, col = qq_col[i],pch=1,lwd=3,cex=cex[3])
#print(max(log.Quantiles))
#
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"
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,lwd=3,cex=cex[3])
if(is.null(signal.col)){
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = t(col)[i],pch=signal.pch[1],cex=signal.cex[1])
}else{
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = signal.col[1],pch=signal.pch[1],cex=signal.cex[1])
}
}
}
}
}
}
box()
if(file.output) grDevices::dev.off()
}
}else{
for(i in 1:R){
print(paste("Q_Q Plotting ",taxa[i],"...",sep=""))
if(file.output){
if(file=="jpg") grDevices::jpeg(paste("QQplot.",taxa[i],".jpg",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi,quality = 100)
if(file=="pdf") grDevices::pdf(paste("GAPIT.Association.Q_Q.",taxa[i],".pdf",sep=""), width = 5.5,height=5.5)
if(file=="tiff") grDevices::tiff(paste("QQplot.",taxa[i],".tiff",sep=""), width = 5.5*dpi,height=5.5*dpi,res=dpi)
graphics::par(mar = c(5,5,4,2),xpd=TRUE)
}else{
if(is.null(grDevices::dev.list())) grDevices::dev.new(width = 5.5, height = 5.5)
graphics::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] <- stats::qbeta(0.95,xi,N-xi+1)
c05[j] <- stats::qbeta(0.05,xi,N-xi+1)
}
index=length(c95):1
}else{
c05 <- 1
c95 <- 1
}
#print(max(log.Quantiles))
#print("@@@@@")
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]))
graphics::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)
graphics::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
#print(log.Quantiles[index])
qq_col = grDevices::rainbow(R)
#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(conf.int) graphics::polygon(c(log.Quantiles[index],log.Quantiles),c(-log10(c05)[index],-log10(c95)),col=qq_col[i],border=conf.int.col)
if(!is.null(threshold.col)){
graphics::par(xpd=FALSE);
graphics::abline(a = 0, b = 1, col = threshold.col[1],lwd=2);
graphics::par(xpd=TRUE)
}
graphics::points(log.Quantiles, log.P.values, col = col[1],pch=19,cex=2)
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){
#print("!!!!")
#cover the points that exceed the threshold with the color "white"
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index], col = "white",pch=19,lwd=3,cex=cex[3])
if(is.null(signal.col)){
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = col[1],pch=signal.pch[1],cex=signal.cex[1])
}else{
graphics::points(log.Quantiles[thre.index],log.P.values[thre.index],col = signal.col[1],pch=signal.pch[1],cex=signal.cex[1])
}
}
}
}
}
box()
if(file.output) grDevices::dev.off()
}
}
print("Multiple QQ plot has been finished!")
}
}#End of Whole function
#}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.