X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FDNA.R;h=722f7b4b9c9cabe100f3a7cd30deed790ff354cb;hb=a8d86ed3981c64b6d20316d54d3ed4d02c2b887e;hp=919f0ecd88912e69ff85b2ee04181d122936eb7b;hpb=8883719a9139685f26d4c0c4cb26872d5f6d0d96;p=ape.git diff --git a/R/DNA.R b/R/DNA.R index 919f0ec..722f7b4 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -1,12 +1,19 @@ -## DNA.R (2009-09-06) +## DNA.R (2010-10-19) ## Manipulations and Comparisons of DNA Sequences -## Copyright 2002-2009 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) { @@ -49,7 +56,7 @@ 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 @@ -84,12 +91,29 @@ 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(...) 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.") @@ -106,6 +130,9 @@ cbind.DNAbin <- 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)) for (i in 1:n) class(obj[[i]]) <- NULL if (check.names) { @@ -144,29 +171,24 @@ cbind.DNAbin <- } c.DNAbin <- function(..., recursive = FALSE) - structure(NextMethod("c"), class = "DNAbin") - -print.DNAbin <- function(x, ...) { - 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) @@ -183,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" @@ -196,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") @@ -269,6 +291,37 @@ base.freq <- function(x, freq = FALSE) 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) @@ -286,33 +339,21 @@ seg.sites <- function(x) 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", "N") - 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 } @@ -321,7 +362,9 @@ 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, n, s, rep(1L, s), PACKAGE = "ape")[[4]]