bug fix in root()

[ape.git] / R / DNA.R

diff --git a/R/DNA.R b/R/DNA.R
index bd5f1954492527573a7aee1780898aef9584904b..e51f252eb7b77421fad835d13ae7671ddae32d00 100644 (file)
--- a/R/DNA.R
+++ b/R/DNA.R
@@ -1,12 +1,33 @@
-## DNA.R (2007-12-21)
+## DNA.R (2008-06-08)
 
 ##   Manipulations and Comparisons of DNA Sequences
 
-## Copyright 2002-2007 Emmanuel Paradis
+## Copyright 2002-2008 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)
+        x[-i]
+    }
+
+    if (class(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, ]
+        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)
@@ -209,25 +230,22 @@ base.freq <- function(x)
 {
     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),
+             DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
     names(BF) <- letters[c(1, 3, 7, 20)]
     BF
 }
 
-GC.content <- function(x)
-{
-    BF <- base.freq(x)
-    sum(BF[2:3])
-}
+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]
-    ans <- .C("SegSites", x, as.integer(n), as.integer(s),
-              integer(s), PACKAGE = "ape")
+    ans <- .C("SegSites", x, n, s, integer(s),
+              DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
     which(as.logical(ans[[4]]))
 }
 
@@ -235,26 +253,12 @@ 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.")
-
-    n <- dim(x)
-    s <- n[2]
-    n <- n[1]
-
-    ## <FIXME> this should be safely deleted
-    if (!pairwise.deletion) {
-        keep <- .C("GlobalDeletionDNA", x, as.integer(n),
-                   as.integer(s), as.integer(rep(1, s)),
-                   PACKAGE = "ape")[[4]]
-        x <- x[,  as.logical(keep)]
-        s <- dim(x)[2]
-    }
-    ## </FIXME>
-
-    ans <- .C("NucleotideDiversity", x, as.integer(n), as.integer(s),
-              as.integer(pairwise.deletion), double(1), PACKAGE = "ape")[[5]]
-
+    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))
+        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
@@ -282,9 +286,8 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
     n <- n[1]
     BF <- if (is.null(base.freq)) base.freq(x) else base.freq
     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]
     }
@@ -292,10 +295,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) {