X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FDNA.R;h=722f7b4b9c9cabe100f3a7cd30deed790ff354cb;hb=2741f6e9f61e33c7b499f27c47604606d08f4bea;hp=b5f472be1ed6f53115220335dee71e00747f2052;hpb=6dfbab243973c0c3fa2e6d02b190aefbe5a67280;p=ape.git

diff --git a/R/DNA.R b/R/DNA.R
index b5f472b..722f7b4 100644
--- a/R/DNA.R
+++ b/R/DNA.R
@@ -1,12 +1,41 @@
-## DNA.R (2008-01-19)
+## DNA.R (2010-10-19)
 
 ##   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)
+        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)
@@ -27,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]
@@ -42,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.")
@@ -61,67 +91,113 @@ as.matrix.DNAbin <- function(x, ...)
     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(...)
-    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]
@@ -129,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"
@@ -142,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")
@@ -205,45 +281,62 @@ 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", as.raw(x), as.integer(n),
-             double(4), PACKAGE = "ape")[[3]]
+    BF <- .C("BaseProportion", x, n, double(4), freq,
+             DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
     names(BF) <- letters[c(1, 3, 7, 20)]
     BF
 }
 
-GC.content <- function(x)
+Ftab <- function(x, y = NULL)
 {
-    BF <- base.freq(x)
-    sum(BF[2:3])
+    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
 }
 
-seg.sites <- function(x)
-{
-    n <- dim(x)
-    s <- n[2]
-    n <- n[1]
-    ans <- .C("SegSites", x, as.integer(n), as.integer(s),
-              integer(s), PACKAGE = "ape")
-    which(as.logical(ans[[4]]))
-}
+GC.content <- function(x) sum(base.freq(x)[2:3])
 
-nuc.div <- function(x, variance = FALSE, pairwise.deletion = FALSE)
+seg.sites <- function(x)
 {
-    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)*s) + 2*(n^2 + n + 3)*ans/(9*n*(n - 1))
-        ans <- c(ans, var)
+    if (is.vector(x)) n <- 1
+    else { # 'x' is a matrix
+        n <- dim(x)
+        s <- n[2]
+        n <- n[1]
     }
-    ans
+    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]]))
 }
 
 dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
@@ -251,13 +344,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")
-    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.")
         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
     }
@@ -266,11 +362,12 @@ 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, as.integer(n),
-                   as.integer(s), as.integer(rep(1, s)),
-                   PACKAGE = "ape")[[4]]
+        keep <- .C("GlobalDeletionDNA", x, n, s,
+                   rep(1L, s), PACKAGE = "ape")[[4]]
         x <- x[,  as.logical(keep)]
         s <- dim(x)[2]
     }
@@ -278,10 +375,10 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
     var <- if (variance) double(Ndist) else 0
     if (!gamma) gamma <- alpha <- 0
     else alpha <- gamma <- 1
-    d <- .C("dist_dna", x, as.integer(n), as.integer(s),
-            as.integer(imod), double(Ndist), BF,
+    d <- .C("dist_dna", x, n, s, imod, double(Ndist), BF,
             as.integer(pairwise.deletion), as.integer(variance),
-            var, as.integer(gamma), alpha, PACKAGE = "ape")
+            var, as.integer(gamma), alpha, DUP = FALSE, NAOK = TRUE,
+            PACKAGE = "ape")
     if (variance) var <- d[[9]]
     d <- d[[5]]
     if (imod == 11) {