X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FDNA.R;h=48e605484f877964815eef05fe055bbe51a69587;hb=a0436318d70829a2d16134be7ca1d6d454613a20;hp=90e841badb278b4aed61031d6df4f4a50cfe14d2;hpb=41cdf69dd7d0b334a94575c857f0d4c91a61e86f;p=ape.git diff --git a/R/DNA.R b/R/DNA.R index 90e841b..48e6054 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -1,24 +1,32 @@ -## DNA.R (2008-08-07) +## DNA.R (2013-01-04) ## Manipulations and Comparisons of DNA Sequences -## Copyright 2002-2008 Emmanuel Paradis +## Copyright 2002-2013 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,102 +56,135 @@ as.alignment <- function(x) obj } -"[.DNAbin" <- function(x, i, j, drop = TRUE) +"[.DNAbin" <- function(x, i, j, drop = FALSE) { - class(x) <- NULL - if (is.matrix(x)) { - if (nargs() == 2 && !missing(i)) ans <- x[i] - else { - nd <- dim(x) - if (missing(i)) i <- 1:nd[1] - if (missing(j)) j <- 1:nd[2] - ans <- x[i, j, drop = drop] - } - } else { - if (missing(i)) i <- 1:length(x) - ans <- x[i] - } - structure(ans, class = "DNAbin") + ans <- NextMethod("[", drop = drop) + class(ans) <- "DNAbin" + 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.") - nms <- names(x) - n <- length(x) - s <- length(x[[1]]) - x <- matrix(unlist(x), n, s, byrow = TRUE) - rownames(x) <- nms - class(x) <- "DNAbin" + if (is.vector(x)) { + dim(x) <- c(1, length(x)) + return(x) } - x + s <- unique(unlist(lapply(x, length))) + if (length(s) != 1) + stop("DNA sequences in list not of the same length.") + nms <- names(x) + n <- length(x) + y <- matrix(as.raw(0), n, s) + for (i in seq_len(n)) y[i, ] <- x[[i]] + rownames(y) <- nms + class(y) <- "DNAbin" + y +} + +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") + nTot <- length(x[[1]]) + cat("Sequence length:", nTot, "\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)) - mini <- min(tmp) - maxi <- max(tmp) + tmp <- unlist(lapply(x, length)) + nTot <- sum(tmp) + mini <- range(tmp) + maxi <- mini[2] + mini <- mini[1] if (mini == maxi) cat("All sequences of same length:", maxi, "\n") else { @@ -158,30 +199,37 @@ 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) - 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" - if (printlen < nd[1]) { - nms <- nms[1:printlen] - TAIL <- "...\n\n" - } - 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") + nTot <- length(x) + 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" + if (printlen < nd[1]) { + nms <- nms[1:printlen] + TAIL <- "...\n\n" + } + cat("\nLabels:", paste(nms, collapse = " "), TAIL) + } else { + cat("1 DNA sequence in binary format stored in a vector.\n\n") + cat("Sequence length:", nTot, "\n\n") + } } - cat("Base composition:\n") - print(round(base.freq(object), digits)) + if (nTot <= 1e6) { + cat("Base composition:\n") + print(round(base.freq(x), digits)) + } else cat("More than 1 million nucleotides: not printing base composition\n") } as.DNAbin <- function(x, ...) UseMethod("as.DNAbin") -._cs_<- letters[c(1, 7, 3, 20, 18, 13, 23, 19, 11, 25, 22, 8, 4, 2, 14)] +._cs_ <- c("a", "g", "c", "t", "r", "m", "w", "s", "k", + "y", "v", "h", "d", "b", "n", "-", "?") -._bs_<- c(136, 72, 40, 24, 192, 160, 144, 96, 80, 48, 224, 176, 208, 112, 240) +._bs_ <- c(136, 72, 40, 24, 192, 160, 144, 96, 80, + 48, 224, 176, 208, 112, 240, 4, 2) as.DNAbin.character <- function(x, ...) { @@ -234,42 +282,78 @@ 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, all = FALSE) { - if (is.list(x)) x <- unlist(x) - n <- length(x) - BF <- .C("BaseProportion", x, n, double(4), - DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]] - names(BF) <- letters[c(1, 3, 7, 20)] + f <- function(x) + .C("BaseProportion", x, as.integer(length(x)), double(17), + DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]] + + if (is.list(x)) { + BF <- rowSums(sapply(x, f)) + n <- sum(sapply(x, length)) + } else { + n <- length(x) + BF <-.C("BaseProportion", x, as.integer(n), double(17), + DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]] + } + + names(BF) <- c("a", "c", "g", "t", "r", "m", "w", "s", + "k", "y", "v", "h", "d", "b", "n", "-", "?") + if (all) { + if (!freq) BF <- BF / n + } else { + BF <- BF[1:4] + if (!freq) BF <- BF / sum(BF) + } BF } -GC.content <- function(x) sum(base.freq(x)[2:3]) - -seg.sites <- function(x) +Ftab <- function(x, y = NULL) { - if (is.list(x)) x <- as.matrix(x) - n <- dim(x) - s <- n[2] - n <- n[1] - ans <- .C("SegSites", x, n, s, integer(s), - DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") - which(as.logical(ans[[4]])) + 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 } -nuc.div <- function(x, variance = FALSE, pairwise.deletion = FALSE) +GC.content <- function(x) sum(base.freq(x)[2:3]) + +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)*dim(x)[2]) + 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, as.integer(n), as.integer(s), + integer(s), DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") + which(as.logical(ans[[4]])) } dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, @@ -277,13 +361,17 @@ 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", + "INDEL", "INDELBLOCK") + 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.") + warning("computing variance 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:17)) { warning(paste("gamma-correction not available for model", model)) gamma <- FALSE } @@ -292,21 +380,30 @@ 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 (imod %in% 16:17) pairwise.deletion <- TRUE + if (!pairwise.deletion) { keep <- .C("GlobalDeletionDNA", x, n, s, rep(1L, s), PACKAGE = "ape")[[4]] - x <- x[, as.logical(keep)] + x <- x[, as.logical(keep)] s <- dim(x)[2] } Ndist <- if (imod == 11) n*n else n*(n - 1)/2 var <- if (variance) double(Ndist) else 0 if (!gamma) gamma <- alpha <- 0 - else alpha <- gamma <- 1 - d <- .C("dist_dna", x, n, s, imod, double(Ndist), BF, - as.integer(pairwise.deletion), as.integer(variance), - var, as.integer(gamma), alpha, DUP = FALSE, NAOK = TRUE, - PACKAGE = "ape") + else { + alpha <- gamma + gamma <- 1 + } + d <- .C("dist_dna", x, as.integer(n), as.integer(s), imod, + double(Ndist), BF, as.integer(pairwise.deletion), + as.integer(variance), var, as.integer(gamma), + alpha, DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") if (variance) var <- d[[9]] d <- d[[5]] if (imod == 11) { @@ -324,3 +421,82 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, if (variance) attr(d, "variance") <- var d } + +image.DNAbin <- function(x, what, col, bg = "white", xlab = "", ylab = "", + show.labels = TRUE, cex.lab = 1, legend = TRUE, ...) +{ + what <- + if (missing(what)) c("a", "g", "c", "t", "n", "-") else tolower(what) + if (missing(col)) + col <- c("red", "yellow", "green", "blue", "grey", "black") + n <- (dx <- dim(x))[1] # number of sequences + s <- dx[2] # number of sites + y <- integer(N <- length(x)) + ncl <- length(what) + col <- rep(col, length.out = ncl) + brks <- 0.5:(ncl + 0.5) + sm <- 0L + for (i in ncl:1) { + k <- ._bs_[._cs_ == what[i]] + sel <- which(x == k) + if (L <- length(sel)) { + y[sel] <- i + sm <- sm + L + } else { + what <- what[-i] + col <- col[-i] + brks <- brks[-i] + } + } + dim(y) <- dx + ## if there's no 0 in y, must drop 'bg' from the cols passed to image: + if (sm == N) { + leg.co <- co <- col + leg.txt <- toupper(what) + } else { + co <- c(bg, col) + leg.txt <- c(toupper(what), "others") + leg.co <- c(col, bg) + brks <- c(-0.5, brks) + } + yaxt <- if (show.labels) "n" else "s" + graphics::image.default(1:s, 1:n, t(y), col = co, xlab = xlab, + ylab = ylab, yaxt = yaxt, breaks = brks, ...) + if (show.labels) + mtext(rownames(x), side = 2, line = 0.1, at = 1:n, + cex = cex.lab, adj = 1, las = 1) + if (legend) { + psr <- par("usr") + xx <- psr[2]/2 + yy <- psr[4] * (0.5 + 0.5/par("plt")[4]) + legend(xx, yy, legend = leg.txt, pch = 22, pt.bg = leg.co, + pt.cex = 2, bty = "n", xjust = 0.5, yjust = 0.5, + horiz = TRUE, xpd = TRUE) + } +} + +where <- function(x, pattern) +{ + pat <- as.DNAbin(strsplit(pattern, NULL)[[1]]) + p <- as.integer(length(pat)) + + foo <- function(x, pat, p) { + s <- as.integer(length(x)) + if (s < p) stop("sequence shorter than the pattern") + ans <- .C("where", x, pat, s, p, integer(s), 0L, + DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") + n <- ans[[6]] + if (n) ans[[5]][seq_len(n)] - p + 2L else integer() + } + + if (is.list(x)) return(lapply(x, foo, pat = pat, p = p)) + if (is.matrix(x)) { + n <- nrow(x) + res <- vector("list", n) + for (i in seq_along(n)) + res[[i]] <- foo(x[i, , drop = TRUE], pat, p) + names(res) <- rownames(x) + return(res) + } + foo(x, pat, p) # if x is a vector +}