3 ## Manipulations and Comparisons of DNA Sequences
5 ## Copyright 2002-2013 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 labels.DNAbin <- function(object, ...)
12 if (is.list(object)) return(names(object))
13 if (is.matrix(object)) return(rownames(object))
17 del.gaps <- function(x)
19 deleteGaps <- function(x) {
21 if (length(i)) x[-i] else x
24 if (!inherits(x, "DNAbin")) x <- as.DNAbin(x)
27 y <- vector("list", n)
28 for (i in 1:n) y[[i]] <- x[i, ]
29 names(y) <- rownames(x)
33 if (!is.list(x)) return(deleteGaps(x))
34 x <- lapply(x, deleteGaps)
39 as.alignment <- function(x)
41 if (is.list(x)) n <- length(x)
42 if (is.matrix(x)) n <- dim(x)[1]
47 seq[i] <- paste(x[[i]], collapse = "")
50 nam <- dimnames(x)[[1]]
52 seq[i] <- paste(x[i, ], collapse = "")
54 obj <- list(nb = n, seq = seq, nam = nam, com = NA)
55 class(obj) <- "alignment"
59 "[.DNAbin" <- function(x, i, j, drop = FALSE)
61 ans <- NextMethod("[", drop = drop)
62 class(ans) <- "DNAbin"
66 as.matrix.DNAbin <- function(x, ...)
68 if (is.matrix(x)) return(x)
70 dim(x) <- c(1, length(x))
73 s <- unique(unlist(lapply(x, length)))
75 stop("DNA sequences in list not of the same length.")
78 y <- matrix(as.raw(0), n, s)
79 for (i in seq_len(n)) y[i, ] <- x[[i]]
85 as.list.DNAbin <- function(x, ...)
87 if (is.list(x)) return(x)
88 if (is.null(dim(x))) obj <- list(x) # cause is.vector() doesn't work
91 obj <- vector("list", n)
92 for (i in 1:n) obj[[i]] <- x[i, ]
93 names(obj) <- rownames(x)
95 class(obj) <- "DNAbin"
99 rbind.DNAbin <- function(...)
100 ### works only with matrices for the moment
104 if (n == 1) return(obj[[1]])
106 if (!is.matrix(obj[[1]]))
107 stop("the 'rbind' method for \"DNAbin\" accepts only matrices")
108 NC <- unlist(lapply(obj, ncol))
109 if (length(unique(NC)) > 1)
110 stop("matrices do not have the same number of columns.")
111 for (i in 1:n) class(obj[[i]]) <- NULL
112 for (i in 2:n) obj[[1]] <- rbind(obj[[1]], obj[[i]])
113 structure(obj[[1]], class = "DNAbin")
117 function(..., check.names = TRUE, fill.with.gaps = FALSE,
119 ### works only with matrices for the moment
123 if (n == 1) return(obj[[1]])
125 if (!is.matrix(obj[[1]]))
126 stop("the 'cbind' method for \"DNAbin\" accepts only matrices")
127 NR <- unlist(lapply(obj, nrow))
128 for (i in 1:n) class(obj[[i]]) <- NULL
130 nms <- unlist(lapply(obj, rownames))
131 if (fill.with.gaps) {
132 NC <- unlist(lapply(obj, ncol))
134 ans <- matrix(as.raw(4), length(nms), sum(NC))
138 to <- from + NC[i] - 1
139 tmp <- rownames(obj[[i]])
140 nmsi <- tmp[tmp %in% nms]
141 ans[nmsi, from:to] <- obj[[i]][nmsi, , drop = FALSE]
147 nms <- names(tab)[which(ubi)]
148 ans <- obj[[1]][nms, , drop = FALSE]
150 ans <- cbind(ans, obj[[i]][nms, , drop = FALSE])
151 if (!quiet && !all(ubi))
152 warning("some rows were dropped.")
155 if (length(unique(NR)) > 1)
156 stop("matrices do not have the same number of rows.")
157 ans <- matrix(unlist(obj), NR)
158 rownames(ans) <- rownames(obj[[1]])
160 class(ans) <- "DNAbin"
164 c.DNAbin <- function(..., recursive = FALSE)
166 if (!all(unlist(lapply(list(...), is.list))))
167 stop("the 'c' method for \"DNAbin\" accepts only lists")
168 structure(NextMethod("c"), class = "DNAbin")
171 print.DNAbin <- function(x, printlen = 6, digits = 3, ...)
177 cat("1 DNA sequence in binary format stored in a list.\n\n")
178 nTot <- length(x[[1]])
179 cat("Sequence length:", nTot, "\n\n")
180 cat("Label:", nms, "\n\n")
182 cat(n, "DNA sequences in binary format stored in a list.\n\n")
183 tmp <- unlist(lapply(x, length))
189 cat("All sequences of same length:", maxi, "\n")
191 cat("Mean sequence length:", round(mean(tmp), 3), "\n")
192 cat(" Shortest sequence:", mini, "\n")
193 cat(" Longest sequence:", maxi, "\n")
197 nms <- nms[1:printlen]
200 cat("\nLabels:", paste(nms, collapse = " "), TAIL)
207 cat(nd[1], "DNA sequences in binary format stored in a matrix.\n\n")
208 cat("All sequences of same length:", nd[2], "\n")
210 if (printlen < nd[1]) {
211 nms <- nms[1:printlen]
214 cat("\nLabels:", paste(nms, collapse = " "), TAIL)
216 cat("1 DNA sequence in binary format stored in a vector.\n\n")
217 cat("Sequence length:", nTot, "\n\n")
221 cat("Base composition:\n")
222 print(round(base.freq(x), digits))
223 } else cat("More than 1 million nucleotides: not printing base composition\n")
226 as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
228 ._cs_ <- c("a", "g", "c", "t", "r", "m", "w", "s", "k",
229 "y", "v", "h", "d", "b", "n", "-", "?")
231 ._bs_ <- c(136, 72, 40, 24, 192, 160, 144, 96, 80,
232 48, 224, 176, 208, 112, 240, 4, 2)
234 as.DNAbin.character <- function(x, ...)
239 ans[which(x == ._cs_[i])] <- as.raw(._bs_[i])
240 ans[which(x == "-")] <- as.raw(4)
241 ans[which(x == "?")] <- as.raw(2)
244 dimnames(ans) <- dimnames(x)
246 class(ans) <- "DNAbin"
250 as.DNAbin.alignment <- function(x, ...)
253 x$seq <- tolower(x$seq)
254 ans <- matrix("", n, nchar(x$seq[1]))
256 ans[i, ] <- strsplit(x$seq[i], "")[[1]]
257 rownames(ans) <- gsub(" +$", "", gsub("^ +", "", x$nam))
258 as.DNAbin.character(ans)
261 as.DNAbin.list <- function(x, ...)
263 obj <- lapply(x, as.DNAbin)
264 class(obj) <- "DNAbin"
268 as.character.DNAbin <- function(x, ...)
271 ans <- character(length(xx))
273 ans[which(xx == ._bs_[i])] <- ._cs_[i]
274 ans[which(xx == 4)] <- "-"
275 ans[which(xx == 2)] <- "?"
278 dimnames(ans) <- dimnames(xx)
282 if (is.list(x)) lapply(x, f) else f(x)
285 base.freq <- function(x, freq = FALSE, all = FALSE)
288 .C("BaseProportion", x, as.integer(length(x)), double(17),
289 DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
292 BF <- rowSums(sapply(x, f))
293 n <- sum(sapply(x, length))
296 BF <-.C("BaseProportion", x, as.integer(n), double(17),
297 DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
300 names(BF) <- c("a", "c", "g", "t", "r", "m", "w", "s",
301 "k", "y", "v", "h", "d", "b", "n", "-", "?")
303 if (!freq) BF <- BF / n
306 if (!freq) BF <- BF / sum(BF)
311 Ftab <- function(x, y = NULL)
317 if (length(x) != length(y))
318 stop("'x' and 'y' not of same lenght")
319 } else { # 'x' is a matrix
320 y <- x[2, , drop = TRUE]
321 x <- x[1, , drop = TRUE]
326 if (length(x) != length(y))
327 stop("'x' and 'y' not of same lenght")
329 out <- matrix(0, 4, 4)
330 k <- c(136, 40, 72, 24)
335 out[i, j] <- sum(a & b)
338 dimnames(out)[1:2] <- list(c("a", "c", "g", "t"))
342 GC.content <- function(x) sum(base.freq(x)[2:3])
344 seg.sites <- function(x)
346 if (is.list(x)) x <- as.matrix(x)
347 if (is.vector(x)) n <- 1
348 else { # 'x' is a matrix
353 if (n == 1) return(integer(0))
354 ans <- .C("SegSites", x, as.integer(n), as.integer(s),
355 integer(s), DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
356 which(as.logical(ans[[4]]))
359 dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
360 pairwise.deletion = FALSE, base.freq = NULL,
363 MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93",
364 "GG95", "LOGDET", "BH87", "PARALIN", "N", "TS", "TV",
365 "INDEL", "INDELBLOCK")
366 imod <- pmatch(toupper(model), MODELS)
368 stop(paste("'model' must be one of:",
369 paste("\"", MODELS, "\"", sep = "", collapse = " ")))
370 if (imod == 11 && variance) {
371 warning("computing variance not available for model BH87")
374 if (gamma && imod %in% c(1, 5:7, 9:17)) {
375 warning(paste("gamma-correction not available for model", model))
378 if (is.list(x)) x <- as.matrix(x)
379 nms <- dimnames(x)[[1]]
384 if (imod %in% c(4, 6:8)) {
385 BF <- if (is.null(base.freq)) base.freq(x) else base.freq
388 if (imod %in% 16:17) pairwise.deletion <- TRUE
390 if (!pairwise.deletion) {
391 keep <- .C("GlobalDeletionDNA", x, n, s,
392 rep(1L, s), PACKAGE = "ape")[[4]]
393 x <- x[, as.logical(keep)]
396 Ndist <- if (imod == 11) n*n else n*(n - 1)/2
397 var <- if (variance) double(Ndist) else 0
398 if (!gamma) gamma <- alpha <- 0
403 d <- .C("dist_dna", x, as.integer(n), as.integer(s), imod,
404 double(Ndist), BF, as.integer(pairwise.deletion),
405 as.integer(variance), var, as.integer(gamma),
406 alpha, DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
407 if (variance) var <- d[[9]]
411 dimnames(d) <- list(nms, nms)
414 attr(d, "Labels") <- nms
415 attr(d, "Diag") <- attr(d, "Upper") <- FALSE
416 attr(d, "call") <- match.call()
417 attr(d, "method") <- model
419 if (as.matrix) d <- as.matrix(d)
421 if (variance) attr(d, "variance") <- var
425 image.DNAbin <- function(x, what, col, bg = "white", xlab = "", ylab = "",
426 show.labels = TRUE, cex.lab = 1, legend = TRUE, ...)
429 if (missing(what)) c("a", "g", "c", "t", "n", "-") else tolower(what)
431 col <- c("red", "yellow", "green", "blue", "grey", "black")
432 n <- (dx <- dim(x))[1] # number of sequences
433 s <- dx[2] # number of sites
434 y <- integer(N <- length(x))
436 col <- rep(col, length.out = ncl)
437 brks <- 0.5:(ncl + 0.5)
440 k <- ._bs_[._cs_ == what[i]]
442 if (L <- length(sel)) {
452 ## if there's no 0 in y, must drop 'bg' from the cols passed to image:
455 leg.txt <- toupper(what)
458 leg.txt <- c(toupper(what), "others")
460 brks <- c(-0.5, brks)
462 yaxt <- if (show.labels) "n" else "s"
463 graphics::image.default(1:s, 1:n, t(y), col = co, xlab = xlab,
464 ylab = ylab, yaxt = yaxt, breaks = brks, ...)
466 mtext(rownames(x), side = 2, line = 0.1, at = 1:n,
467 cex = cex.lab, adj = 1, las = 1)
471 yy <- psr[4] * (0.5 + 0.5/par("plt")[4])
472 legend(xx, yy, legend = leg.txt, pch = 22, pt.bg = leg.co,
473 pt.cex = 2, bty = "n", xjust = 0.5, yjust = 0.5,
474 horiz = TRUE, xpd = TRUE)
478 where <- function(x, pattern)
480 pat <- as.DNAbin(strsplit(pattern, NULL)[[1]])
481 p <- as.integer(length(pat))
483 foo <- function(x, pat, p) {
484 s <- as.integer(length(x))
485 if (s < p) stop("sequence shorter than the pattern")
486 ans <- .C("where", x, pat, s, p, integer(s), 0L,
487 DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
489 if (n) ans[[5]][seq_len(n)] - p + 2L else integer()
492 if (is.list(x)) return(lapply(x, foo, pat = pat, p = p))
495 res <- vector("list", n)
496 for (i in seq_along(n))
497 res[[i]] <- foo(x[i, , drop = TRUE], pat, p)
498 names(res) <- rownames(x)
501 foo(x, pat, p) # if x is a vector