-## DNA.R (2011-02-18)
+## DNA.R (2011-03-21)
## Manipulations and Comparisons of DNA Sequences
as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
-._cs_<- letters[c(1, 7, 3, 20, 18, 13, 23, 19, 11, 25, 22, 8, 4, 2, 14)]
+._cs_ <- c("a", "g", "c", "t", "r", "m", "w", "s", "k",
+ "y", "v", "h", "d", "b", "n", "-", "?")
-._bs_<- c(136, 72, 40, 24, 192, 160, 144, 96, 80, 48, 224, 176, 208, 112, 240)
+._bs_ <- c(136, 72, 40, 24, 192, 160, 144, 96, 80,
+ 48, 224, 176, 208, 112, 240, 4, 2)
as.DNAbin.character <- function(x, ...)
{
stop(paste("'model' must be one of:",
paste("\"", MODELS, "\"", sep = "", collapse = " ")))
if (imod == 11 && variance) {
- warning("computing variance temporarily not available for model BH87.")
+ warning("computing variance not available for model BH87.")
variance <- FALSE
}
if (gamma && imod %in% c(1, 5:7, 9:15)) {
if (variance) attr(d, "variance") <- var
d
}
+
+image.DNAbin <- function(x, what, col, bg = "white", xlab = "", ylab = "",
+ show.labels = TRUE, cex.lab = 1, legend = TRUE, ...)
+{
+ what <-
+ if (missing(what)) c("a", "g", "c", "t", "n", "-") else tolower(what)
+ if (missing(col))
+ col <- c("red", "yellow", "green", "blue", "grey", "black")
+ n <- (dx <- dim(x))[1] # number of sequences
+ s <- dx[2] # number of sites
+ y <- integer(N <- length(x))
+ ncl <- length(what)
+ col <- rep(col, length.out = ncl)
+ brks <- 0.5:(ncl + 0.5)
+ sm <- 0L
+ for (i in ncl:1) {
+ k <- ._bs_[._cs_ == what[i]]
+ sel <- which(x == k)
+ if (L <- length(sel)) {
+ y[sel] <- i
+ sm <- sm + L
+ } else {
+ what <- what[-i]
+ col <- col[-i]
+ brks <- brks[-i]
+ }
+ }
+ dim(y) <- dx
+ ## if there's no 0 in y, must drop 'bg' from the cols passed to image:
+ if (sm == N) {
+ leg.co <- co <- col
+ leg.txt <- toupper(what)
+ } else {
+ co <- c(bg, col)
+ leg.txt <- c(toupper(what), "others")
+ leg.co <- c(col, bg)
+ brks <- c(-0.5, brks)
+ }
+ yaxt <- if (show.labels) "n" else "s"
+ graphics::image.default(1:s, 1:n, t(y), col = co, xlab = xlab,
+ ylab = ylab, yaxt = yaxt, breaks = brks, ...)
+ if (show.labels)
+ mtext(rownames(x), side = 2, line = 0.1, at = 1:n,
+ cex = cex.lab, adj = 1, las = 1)
+ if (legend) {
+ psr <- par("usr")
+ xx <- psr[2]/2
+ yy <- psr[4] * (0.5 + 0.5/par("plt")[4])
+ legend(xx, yy, legend = leg.txt, pch = 22, pt.bg = leg.co,
+ pt.cex = 2, bty = "n", xjust = 0.5, yjust = 0.5,
+ horiz = TRUE, xpd = TRUE)
+ }
+}