cleaning of C files and update on simulation of OU process

[ape.git] / R / DNA.R

diff --git a/R/DNA.R b/R/DNA.R
index e71f14224450301526daec4b3621ad61451a5e96..63306ef9cc78e3c2556c8196f71520e931ca625e 100644 (file)
--- a/R/DNA.R
+++ b/R/DNA.R
@@ -1,4 +1,4 @@
-## DNA.R (2011-02-18)
+## DNA.R (2011-03-21)
 
 ##   Manipulations and Comparisons of DNA Sequences
 
@@ -226,9 +226,11 @@ print.DNAbin <- function(x, printlen = 6, digits = 3, ...)
 
 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, ...)
 {
@@ -403,3 +405,56 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
     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)
+    }
+}