sigBarPlot.R – A function that plots significant difference comparator brackets for various levels given their p-values.

sigBarPlot <- function(x,y,sig,spaceFac,yRange,xRangeStar,linetype,color,linewidth){
	# x: x-data, N <- length(x)
	# y: y-data, length(x) must == length(y)
	# sig: a 3-by-M matrix; column 1 is the x-data, column 2 is the corresponding y-data and column 3 in the corresponding p=value; M must be > 0

	sigBars <- function(xi,xf,y11,y21,y2,p,ySpace,xOut){
		segments(x0=xi, y0=y11, x1=xi, y1=y2, lwd=linewidth, lty=linetype, col=color) # add three line segments to make a "comparator" line
		segments(x0=xf, y0=y21, x1=xf, y1=y2, lwd=linewidth, lty=linetype, col=color)
		segments(x0=xi, y0=y2, x1=(xi+xf)/2-xOut, y1=y2, lwd=linewidth, lty=linetype, col=color)
		segments(x0=(xi+xf)/2+xOut, y0=y2, x1=xf, y1=y2, lwd=linewidth, lty=linetype, col=color)
		if (p==1){text(x=(xi+xf)/2, y=y2, labels="*")} else
			{if(p==2){text(x=(xi+xf)/2, y=y2, labels="**")} else text(x=(xi+xf)/2, y=y2, labels="***")}
	}

	x <- as.vector(x)
	xOut <- (max(x) - min(x))/xRangeStar

	#yMax <- max(y)
	#ySpace <- yMax*spaceFac
	ySpace <- spaceFac*yRange

	M <- length(x)
	Ncombs <- length(sig[,1])
	diffY <- mat.or.vec(Ncombs,1)

	for (i in 1:Ncombs){
		if(sig[i,3] < 0.001){sig[i,3] <- 3}
		else if(sig[i,3] < 0.01){sig[i,3] <- 2}
		else if(sig[i,3] < 0.05){sig[i,3] <- 1}
		else{}
		diffY[i] <- abs(y[sig[i,1]]-y[sig[i,2]])
	}

	if (Ncombs > 1){
		index <- sort.int(diffY, partial=NULL, na.last=NA, decreasing=FALSE,
	         method=c("shell"), index.return=TRUE)
	    index <- index$ix
	    indexMin <- index[1]
	    Sig <- sig[indexMin,]
	    sig <- sig[-indexMin,]

		yBase <- y
		y2 <- max(yBase[Sig[1]:Sig[2]]) + 7/4*ySpace

		for (n in 1:Ncombs){

			x1 <-  x[Sig[1]]
			x2 <-  x[Sig[2]]
			y11 <- yBase[Sig[1]] + ySpace
			y21 <- yBase[Sig[2]] + ySpace
			yBase[Sig[1]:Sig[2]] <- y2

			sigBars(x1,x2,y11,y21,y2,Sig[3],ySpace,xOut)

			if ( n < Ncombs ){

			    if (n == Ncombs-1){Sig <- sig}
		    	else{
		    		diffY <- abs(yBase[sig[,1]]-yBase[sig[,2]])
					index <- sort.int(diffY, partial=NULL, na.last=NA, decreasing=FALSE,
				         method=c("shell"), index.return=TRUE)
				    index <- index$ix
					indexMin <- index[1]
					Sig <- sig[indexMin,]
				    sig <- sig[-indexMin,]
		    	}
		    	y2  <- max(c(yBase[Sig[1]:Sig[2]])) + 7/4*ySpace

			}
		}
	}
	else{
		x1 <- x[1]
		x2 <- x[2]
		y11 <- y[1] + ySpace
		y21 <- y[2] + ySpace
		y2  <- max(c(y11,y21)) + 7/4*ySpace
		sigBars(x1,x2,y11,y21,y2,sig[3],ySpace,xOut)
	}
}

Driver script for the above function.

rm(list = ls())
require(gplots)
source( "sigBarPlot.R" )

N <- 6

x <- seq(1,N)
meany <- x + 2*(runif(N)^2)
sey <- runif(N)
y <- meany + sey

sig1 <- c(1,6,0.0006)
sig2 <- c(1,5,0.008)
sig3 <- c(2,6,0.02)

sig <- rbind(sig1,sig2,sig3)

lineType <- "solid"
lineWidth <- 0.75
lineColor <- "gray"

pdf(file="sigPlot.pdf",width=6,height=6,pointsize=10)
par(mar=c(3,3,0,0), oma=c(0,0,0,0))

xaxis <- barplot(meany, beside=TRUE,
	names.arg = as.character(x),
	#space = s,
  	las = 1,
	col = "gray",
	border = "gray",
    cex.names = 1,
    xaxs="r",
    yaxs="r",
    ylim=c(0,1.2*max(y)))
plotCI(x=xaxis, y=meany, uiw=sey, liw=sey, ui=NULL, li=NULL, err='y', lty=1, pch="", sfrac=0.01, gap=0.01, lwd=1, labels=FALSE,  add=TRUE)
sigBarPlotFast(xaxis,y,sig,0.01,1.2*max(y),15,lineType,lineColor,lineWidth)
mtext("x",1,line=2)
mtext("y",2,line=2)

dev.off()

Figure output: