-## DNA.R (2008-10-08)
+## DNA.R (2011-03-21)
## Manipulations and Comparisons of DNA Sequences
-## Copyright 2002-2008 Emmanuel Paradis
+## Copyright 2002-2011 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
+labels.DNAbin <- function(object, ...)
+{
+ if (is.list(object)) return(names(object))
+ if (is.matrix(object)) return(rownames(object))
+ NULL
+}
+
del.gaps <- function(x)
{
deleteGaps <- function(x) {
if (length(i)) x[-i] else x
}
- if (class(x) != "DNAbin") x <- as.DNAbin(x)
+ if (!inherits(x, "DNAbin")) x <- as.DNAbin(x)
if (is.matrix(x)) {
n <- dim(x)[1]
y <- vector("list", n)
for (i in 1:n) y[[i]] <- x[i, ]
+ names(y) <- rownames(x)
x <- y
rm(y)
}
obj
}
-"[.DNAbin" <- function(x, i, j, drop = TRUE)
+"[.DNAbin" <- function(x, i, j, drop = FALSE)
{
+ oc <- oldClass(x)
class(x) <- NULL
if (is.matrix(x)) {
if (nargs() == 2 && !missing(i)) ans <- x[i]
if (missing(i)) i <- 1:length(x)
ans <- x[i]
}
- structure(ans, class = "DNAbin")
+ class(ans) <- oc
+ ans
}
as.matrix.DNAbin <- function(x, ...)
{
- if (is.matrix(x)) return(x)
if (is.list(x)) {
if (length(unique(unlist(lapply(x, length)))) != 1)
stop("DNA sequences in list not of the same length.")
x
}
+as.list.DNAbin <- function(x, ...)
+{
+ if (is.list(x)) return(x)
+ if (is.null(dim(x))) obj <- list(x) # cause is.vector() doesn't work
+ else { # matrix
+ n <- nrow(x)
+ obj <- vector("list", n)
+ for (i in 1:n) obj[[i]] <- x[i, ]
+ names(obj) <- rownames(x)
+ }
+ class(obj) <- "DNAbin"
+ obj
+}
+
rbind.DNAbin <- function(...)
### works only with matrices for the moment
{
obj <- list(...)
n <- length(obj)
if (n == 1) return(obj[[1]])
+ for (i in 1:n)
+ if (!is.matrix(obj[[1]]))
+ stop("the 'rbind' method for \"DNAbin\" accepts only matrices")
NC <- unlist(lapply(obj, ncol))
if (length(unique(NC)) > 1)
stop("matrices do not have the same number of columns.")
structure(obj[[1]], class = "DNAbin")
}
-cbind.DNAbin <- function(..., check.names = TRUE, fill.with.gaps = FALSE,
- quiet = TRUE)
+cbind.DNAbin <-
+ function(..., check.names = TRUE, fill.with.gaps = FALSE,
+ quiet = FALSE)
### works only with matrices for the moment
{
obj <- list(...)
n <- length(obj)
if (n == 1) return(obj[[1]])
+ for (i in 1:n)
+ if (!is.matrix(obj[[1]]))
+ stop("the 'cbind' method for \"DNAbin\" accepts only matrices")
NR <- unlist(lapply(obj, nrow))
- if (length(unique(NR)) > 1)
- stop("matrices do not have the same number of rows.")
for (i in 1:n) class(obj[[i]]) <- NULL
- nms <- rownames(obj[[1]])
if (check.names) {
- for (i in 2:n)
- if (all(rownames(obj[[i]]) %in% nms))
- obj[[i]] <- obj[[i]][nms, ]
- else stop("rownames do not match among matrices.")
+ nms <- unlist(lapply(obj, rownames))
+ if (fill.with.gaps) {
+ NC <- unlist(lapply(obj, ncol))
+ nms <- unique(nms)
+ ans <- matrix(as.raw(4), length(nms), sum(NC))
+ rownames(ans) <- nms
+ from <- 1
+ for (i in 1:n) {
+ to <- from + NC[i] - 1
+ tmp <- rownames(obj[[i]])
+ nmsi <- tmp[tmp %in% nms]
+ ans[nmsi, from:to] <- obj[[i]][nmsi, , drop = FALSE]
+ from <- to + 1
+ }
+ } else {
+ tab <- table(nms)
+ ubi <- tab == n
+ nms <- names(tab)[which(ubi)]
+ ans <- obj[[1]][nms, , drop = FALSE]
+ for (i in 2:n)
+ ans <- cbind(ans, obj[[i]][nms, , drop = FALSE])
+ if (!quiet && !all(ubi))
+ warning("some rows were dropped.")
+ }
+ } else {
+ if (length(unique(NR)) > 1)
+ stop("matrices do not have the same number of rows.")
+ ans <- matrix(unlist(obj), NR)
+ rownames(ans) <- rownames(obj[[1]])
}
- ans <- matrix(unlist(obj), NR)
- rownames(ans) <- nms
- structure(ans, class = "DNAbin")
+ class(ans) <- "DNAbin"
+ ans
}
-print.DNAbin <- function(x, ...)
+c.DNAbin <- function(..., recursive = FALSE)
{
- n <- 1 # <- if is.vector(x)
- if (is.list(x)) n <- length(x)
- else if (is.matrix(x)) n <- dim(x)[1]
- if (n > 1) cat(n, "DNA sequences in binary format.\n")
- else cat("1 DNA sequence in binary format.\n")
+ if (!all(unlist(lapply(list(...), is.list))))
+ stop("the 'c' method for \"DNAbin\" accepts only lists")
+ structure(NextMethod("c"), class = "DNAbin")
}
-summary.DNAbin <- function(object, printlen = 6, digits = 3, ...)
+print.DNAbin <- function(x, printlen = 6, digits = 3, ...)
{
- if (is.list(object)) {
- n <- length(object)
- nms <- names(object)
+ if (is.list(x)) {
+ n <- length(x)
+ nms <- names(x)
if (n == 1) {
cat("1 DNA sequence in binary format stored in a list.\n\n")
- cat("Sequence length:", length(object[[1]]), "\n\n")
+ cat("Sequence length:", length(x[[1]]), "\n\n")
cat("Label:", nms, "\n\n")
} else {
cat(n, "DNA sequences in binary format stored in a list.\n\n")
- tmp <- unlist(lapply(object, length))
+ tmp <- unlist(lapply(x, length))
mini <- min(tmp)
maxi <- max(tmp)
if (mini == maxi)
}
cat("\nLabels:", paste(nms, collapse = " "), TAIL)
}
- } else if (is.matrix(object)) {
- nd <- dim(object)
- nms <- rownames(object)
+ } else if (is.matrix(x)) {
+ nd <- dim(x)
+ nms <- rownames(x)
cat(nd[1], "DNA sequences in binary format stored in a matrix.\n\n")
cat("All sequences of same length:", nd[2], "\n")
TAIL <- "\n\n"
cat("\nLabels:", paste(nms, collapse = " "), TAIL)
} else {
cat("1 DNA sequence in binary format stored in a vector.\n\n")
- cat("Sequence length:", length(object), "\n\n")
+ cat("Sequence length:", length(x), "\n\n")
}
cat("Base composition:\n")
- print(round(base.freq(object), digits))
+ print(round(base.freq(x), digits))
}
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, ...)
{
if (is.list(x)) lapply(x, f) else f(x)
}
-base.freq <- function(x)
+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),
- 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)
{
if (is.list(x)) x <- as.matrix(x)
- n <- dim(x)
- s <- n[2]
- n <- n[1]
+ if (is.vector(x)) n <- 1
+ else { # 'x' is a matrix
+ n <- dim(x)
+ s <- n[2]
+ n <- n[1]
+ }
+ if (n == 1) return(integer(0))
ans <- .C("SegSites", x, n, s, integer(s),
DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
which(as.logical(ans[[4]]))
}
-nuc.div <- function(x, variance = FALSE, pairwise.deletion = FALSE)
-{
- if (pairwise.deletion && variance)
- warning("cannot compute the variance of nucleotidic diversity\nwith pairwise deletion: try 'pairwise.deletion = FALSE' instead.")
- if (is.list(x)) x <- as.matrix(x)
- n <- dim(x)[1]
- ans <- sum(dist.dna(x, "raw", pairwise.deletion = pairwise.deletion))/
- (n*(n - 1)/2)
- if (variance) {
- var <- (n + 1)*ans/(3*(n + 1)*dim(x)[2]) + 2*(n^2 + n + 3)*ans/(9*n*(n - 1))
- ans <- c(ans, var)
- }
- ans
-}
-
dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
pairwise.deletion = FALSE, base.freq = NULL,
as.matrix = FALSE)
{
MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93",
- "GG95", "LOGDET", "BH87", "PARALIN")
- imod <- which(MODELS == toupper(model))
+ "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
}
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]]
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)
+ }
+}
projects / ape.git / blobdiff