[seqfan] Re: Chewing an old bone

[Russ Cox]

This suggestion sounds interesting.  The best way to explore it
would be if someone could implement it in R, and then we can try
both for a little while.

I am including the R program that generates the current graphs below.
The input file read by read.table is the standard b file format,
without any comments or blank lines.

If someone can adapt this program to generate these alternate
plots, then I will be happy to set it up so that we can see both.

Thanks.
Russ


# R CMD BATCH --no-restore --no-save --vanilla -- script.r

infile = "b123.txt"
outfile = "b123.png"
seqname = "Name Of Sequence"

x = read.table(infile)
  
# Remove infinite values
lgl = is.finite(x[,2])
x = cbind(x[,1][lgl], x[,2][lgl])

png(file=outfile, width=600, height=800)
par(mar=c(5,5,2,2)+0.1, mfcol=c(2,1))

k = min(200, nrow(x))
xsub = x[1:k,]
plot (xsub, 
  ylim = c(min(0, min(xsub[,2])), max(xsub[,2])),
  xlab = "n",
  ylab=paste(seqname, "(n)", sep=""),
  bty="l",
  type="n",
  main = paste("Pin plot of", seqname))
segments(xsub[,1], rep(0, k), xsub[,1], xsub[,2])

# Log the data for the scatterplot?

logplot = T
tval = 100.
if (nrow(x) > 10) {
  fv = fivenum(abs(x[,2]))
  if (all(is.finite(fv))) {
    if (((fv[5]-fv[4]) / (fv[2]-fv[1]+1)) < tval) {
      logplot = F
    }
  }
}

if (logplot) {
  v = x[,2]
  if (all((v+1) >= 0)) {
    plot(x[,1], log10(v+1), pch="+", cex=.5,
      xlab = "n",
      ylab=paste("log(", seqname, "(n)+1)", sep=""),
      main=paste("Scatterplot of log(", seqname, "(n)+1)", sep=""),
      bty="l")
  } else {
    plot(x[,1], log10(abs(v)+1), pch="+", cex=.5,
      xlab = "n",
      ylab=paste("log(|", seqname, "(n)|+1)", sep=""),
      main=paste("Scatterplot of log(", seqname, "(n))", sep=""),
      bty="l")
  }
} else {
  plot(x, pch="+", cex=.5,
    xlab = "n", 
    ylab = paste(seqname, "(n)", sep=""),
    main = paste("Scatterplot of ", seqname, "(n)", sep=""),
    bty="l")
}

dev.off()