# [seqfan] Re: Chewing an old bone

Russ Cox rsc at swtch.com
Thu Oct 13 15:29:34 CEST 2011

```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"

# 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-fv) / (fv-fv+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()
```