X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FDNA.R;h=df5b892192977628392406ecf0b52bd9a91ba1cf;hb=d1546ec66ff1a8ea123adefebe14f6316c23705f;hp=8b016caad40dc9acb8014c4889bee52fa38c2bf0;hpb=92ef6853c5cb0e730dad1af9961064239a1d375f;p=ape.git diff --git a/R/DNA.R b/R/DNA.R index 8b016ca..df5b892 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -1,8 +1,8 @@ -## DNA.R (2010-07-13) +## DNA.R (2011-03-21) ## Manipulations and Comparisons of DNA Sequences -## Copyright 2002-2010 Emmanuel Paradis +## Copyright 2002-2011 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -94,7 +94,7 @@ as.matrix.DNAbin <- function(x, ...) as.list.DNAbin <- function(x, ...) { if (is.list(x)) return(x) - if (is.vector(x)) obj <- list(x) + if (is.null(dim(x))) obj <- list(x) # cause is.vector() doesn't work else { # matrix n <- nrow(x) obj <- vector("list", n) @@ -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, ...) { @@ -281,16 +283,54 @@ as.character.DNAbin <- function(x, ...) if (is.list(x)) lapply(x, f) else f(x) } -base.freq <- function(x, freq = FALSE) +base.freq <- function(x, freq = FALSE, all = FALSE) { if (is.list(x)) x <- unlist(x) n <- length(x) - BF <- .C("BaseProportion", x, n, double(4), freq, - DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]] - names(BF) <- letters[c(1, 3, 7, 20)] + BF <-.C("BaseProportion", x, n, double(17), + DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]] + names(BF) <- c("a", "c", "g", "t", "r", "m", "w", "s", + "k", "y", "v", "h", "d", "b", "n", "-", "?") + if (all) { + if (!freq) BF <- BF / n + } else { + BF <- BF[1:4] + if (!freq) BF <- BF / sum(BF) + } BF } +Ftab <- function(x, y = NULL) +{ + if (is.null(y)) { + if (is.list(x)) { + y <- x[[2]] + x <- x[[1]] + if (length(x) != length(y)) + stop("'x' and 'y' not of same lenght") + } else { # 'x' is a matrix + y <- x[2, , drop = TRUE] + x <- x[1, , drop = TRUE] + } + } else { + x <- as.vector(x) + y <- as.vector(y) + if (length(x) != length(y)) + stop("'x' and 'y' not of same lenght") + } + out <- matrix(0, 4, 4) + k <- c(136, 40, 72, 24) + for (i in 1:4) { + a <- x == k[i] + for (j in 1:4) { + b <- y == k[j] + out[i, j] <- sum(a & b) + } + } + dimnames(out)[1:2] <- list(c("a", "c", "g", "t")) + out +} + GC.content <- function(x) sum(base.freq(x)[2:3]) seg.sites <- function(x) @@ -313,16 +353,16 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, as.matrix = FALSE) { MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93", - "GG95", "LOGDET", "BH87", "PARALIN", "N") + "GG95", "LOGDET", "BH87", "PARALIN", "N", "TS", "TV") imod <- pmatch(toupper(model), MODELS) if (is.na(imod)) 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:12)) { + if (gamma && imod %in% c(1, 5:7, 9:15)) { warning(paste("gamma-correction not available for model", model)) gamma <- FALSE } @@ -331,7 +371,9 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, n <- dim(x) s <- n[2] n <- n[1] - BF <- if (is.null(base.freq)) base.freq(x) else base.freq + if (imod %in% c(4, 6:8)) { + BF <- if (is.null(base.freq)) base.freq(x) else base.freq + } else BF <- 0 if (!pairwise.deletion) { keep <- .C("GlobalDeletionDNA", x, n, s, rep(1L, s), PACKAGE = "ape")[[4]] @@ -363,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) + } +}