getwidth <- function(y) { ly <- length(y) ymin <- min(median(y[1:50]), median(y[(ly-49):ly])) ymax <- max(y) xpeak <- round(mean(which(y == ymax))) thresh <- (ymax - ymin) * 0.5 + ymin rwidth <- rle(y[xpeak:ly] >= thresh)$lengths[1] - 1 lwidth <- rle(rev(y[1:xpeak] >= thresh))$lengths[1] - 1 return(c(thresh, x[xpeak - rwidth], x[xpeak + lwidth])) } pdf(commandArgs(trailingOnly=T)[1],height=6,width=6) plot(x,p,type='l',col=c('red'),main=main,xlab='Distance to the middle',ylab='Percentage') lines(x,m,col=c('blue')) legend('topleft',c('forward tags','reverse tags'),lty=c(1,1,1),col=c('red','blue')) r <- getwidth(p) plen <- r[3] - r[2] abline(h=r[1], col="red", lty=2) text(mean(r[2:3]), r[1], r[3] - r[2], col="red") r <- getwidth(m) mlen <- r[3] - r[2] abline(h=r[1], col="blue", lty=2) text(mean(r[2:3]), r[1], r[3] - r[2], col="blue") elen <- round(mean(c(plen, mlen))) text(0, 0, elen) cat(elen) cat('\n')