X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FDNA.R;h=8b016caad40dc9acb8014c4889bee52fa38c2bf0;hb=92ef6853c5cb0e730dad1af9961064239a1d375f;hp=e51f252eb7b77421fad835d13ae7671ddae32d00;hpb=0eb732cf743bc38c934b860fc6df60aa248b3034;p=ape.git diff --git a/R/DNA.R b/R/DNA.R index e51f252..8b016ca 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -1,24 +1,32 @@ -## DNA.R (2008-06-08) +## DNA.R (2010-07-13) ## Manipulations and Comparisons of DNA Sequences -## Copyright 2002-2008 Emmanuel Paradis +## Copyright 2002-2010 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) { i <- which(x == 4) - x[-i] + 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) } @@ -48,8 +56,9 @@ as.alignment <- function(x) 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] @@ -63,12 +72,12 @@ as.alignment <- function(x) 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.") @@ -82,67 +91,113 @@ as.matrix.DNAbin <- function(x, ...) x } +as.list.DNAbin <- function(x, ...) +{ + if (is.list(x)) return(x) + if (is.vector(x)) obj <- list(x) + 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(...) - 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]]) + 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.") - 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]]) + for (i in 1:n) + if (!is.matrix(obj[[1]])) + stop("the 'cbind' method for \"DNAbin\" accepts only matrices") + 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 } -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") - cat("Summary of sequence lengths:\n") - print(summary(unlist(lapply(object, length)))) + tmp <- unlist(lapply(x, 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] @@ -150,9 +205,9 @@ summary.DNAbin <- function(object, printlen = 6, digits = 3, ...) } 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" @@ -163,10 +218,10 @@ summary.DNAbin <- function(object, printlen = 6, digits = 3, ...) 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") @@ -226,11 +281,11 @@ as.character.DNAbin <- function(x, ...) 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 @@ -241,36 +296,28 @@ 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") + 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