-## DNA.R (2008-04-24)
+## DNA.R (2009-10-02)
## Manipulations and Comparisons of DNA Sequences
-## Copyright 2002-2008 Emmanuel Paradis
+## Copyright 2002-2009 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
+del.gaps <- function(x)
+{
+ deleteGaps <- function(x) {
+ i <- which(x == 4)
+ if (length(i)) x[-i] else 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)
+ }
+ if (!is.list(x)) return(deleteGaps(x))
+ x <- lapply(x, deleteGaps)
+ class(x) <- "DNAbin"
+ x
+}
+
as.alignment <- function(x)
{
if (is.list(x)) n <- length(x)
"[.DNAbin" <- function(x, i, j, drop = TRUE)
{
+ 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.")
### works only with matrices for the moment
{
obj <- list(...)
- nobj <- length(obj)
- if (nobj == 1) stop("only one matrix to bind.")
- NC <- ncol(obj[[1]])
- for (i in 2:nobj)
- if(ncol(obj[[i]]) != NC)
+ n <- length(obj)
+ if (n == 1) return(obj[[1]])
+ NC <- unlist(lapply(obj, ncol))
+ if (length(unique(NC)) > 1)
stop("matrices do not have the same number of columns.")
- for (i in 1:nobj) class(obj[[i]]) <- NULL
- ans <- obj[[1]]
- for (i in 2:nobj) ans <- rbind(ans, obj[[i]])
- structure(ans, class = "DNAbin")
+ for (i in 1:n) class(obj[[i]]) <- NULL
+ for (i in 2:n) obj[[1]] <- rbind(obj[[1]], obj[[i]])
+ structure(obj[[1]], class = "DNAbin")
}
-cbind.DNAbin <- function(..., check.names = TRUE)
+cbind.DNAbin <-
+ function(..., check.names = TRUE, fill.with.gaps = FALSE,
+ quiet = FALSE)
### works only with matrices for the moment
{
obj <- list(...)
- nobj <- length(obj)
- if (nobj == 1) stop("only one matrix to bind.")
- NR <- nrow(obj[[1]])
- for (i in 2:nobj)
- if(nrow(obj[[i]]) != NR)
- stop("matrices do not have the same number of rows.")
- for (i in 1:nobj) class(obj[[i]]) <- NULL
- nms <- rownames(obj[[1]])
+ n <- length(obj)
+ if (n == 1) return(obj[[1]])
+ NR <- unlist(lapply(obj, nrow))
+ for (i in 1:n) class(obj[[i]]) <- NULL
if (check.names) {
- for (i in 2:nobj)
- 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
}
+c.DNAbin <- function(..., recursive = FALSE)
+ structure(NextMethod("c"), class = "DNAbin")
+
print.DNAbin <- function(x, ...)
{
n <- 1 # <- if is.vector(x)
cat("Label:", nms, "\n\n")
} else {
cat(n, "DNA sequences in binary format stored in a list.\n\n")
- cat("Summary of sequence lengths:\n")
- print(summary(unlist(lapply(object, length))))
+ tmp <- unlist(lapply(object, length))
+ mini <- min(tmp)
+ maxi <- max(tmp)
+ if (mini == maxi)
+ cat("All sequences of same length:", maxi, "\n")
+ else {
+ cat("Mean sequence length:", round(mean(tmp), 3), "\n")
+ cat(" Shortest sequence:", mini, "\n")
+ cat(" Longest sequence:", maxi, "\n")
+ }
TAIL <- "\n\n"
if (printlen < n) {
nms <- nms[1:printlen]
if (is.list(x)) lapply(x, f) else f(x)
}
-base.freq <- function(x)
+base.freq <- function(x, freq = FALSE)
{
if (is.list(x)) x <- unlist(x)
n <- length(x)
- BF <- .C("BaseProportion", x, n, double(4),
+ BF <- .C("BaseProportion", x, n, double(4), freq,
DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
names(BF) <- letters[c(1, 3, 7, 20)]
BF
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")
+ 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.")
variance <- FALSE
projects / ape.git / blobdiff