3 ## Manipulations and Comparisons of DNA Sequences
5 ## Copyright 2002-2010 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 del.gaps <- function(x)
12 deleteGaps <- function(x) {
14 if (length(i)) x[-i] else x
17 if (!inherits(x, "DNAbin")) x <- as.DNAbin(x)
20 y <- vector("list", n)
21 for (i in 1:n) y[[i]] <- x[i, ]
22 names(y) <- rownames(x)
26 if (!is.list(x)) return(deleteGaps(x))
27 x <- lapply(x, deleteGaps)
32 as.alignment <- function(x)
34 if (is.list(x)) n <- length(x)
35 if (is.matrix(x)) n <- dim(x)[1]
40 seq[i] <- paste(x[[i]], collapse = "")
43 nam <- dimnames(x)[[1]]
45 seq[i] <- paste(x[i, ], collapse = "")
47 obj <- list(nb = n, seq = seq, nam = nam, com = NA)
48 class(obj) <- "alignment"
52 "[.DNAbin" <- function(x, i, j, drop = TRUE)
57 if (nargs() == 2 && !missing(i)) ans <- x[i]
60 if (missing(i)) i <- 1:nd[1]
61 if (missing(j)) j <- 1:nd[2]
62 ans <- x[i, j, drop = drop]
65 if (missing(i)) i <- 1:length(x)
72 as.matrix.DNAbin <- function(x, ...)
75 if (length(unique(unlist(lapply(x, length)))) != 1)
76 stop("DNA sequences in list not of the same length.")
80 x <- matrix(unlist(x), n, s, byrow = TRUE)
87 rbind.DNAbin <- function(...)
88 ### works only with matrices for the moment
92 if (n == 1) return(obj[[1]])
94 if (!is.matrix(obj[[1]]))
95 stop("the 'rbind' method for \"DNAbin\" accepts only matrices")
96 NC <- unlist(lapply(obj, ncol))
97 if (length(unique(NC)) > 1)
98 stop("matrices do not have the same number of columns.")
99 for (i in 1:n) class(obj[[i]]) <- NULL
100 for (i in 2:n) obj[[1]] <- rbind(obj[[1]], obj[[i]])
101 structure(obj[[1]], class = "DNAbin")
105 function(..., check.names = TRUE, fill.with.gaps = FALSE,
107 ### works only with matrices for the moment
111 if (n == 1) return(obj[[1]])
113 if (!is.matrix(obj[[1]]))
114 stop("the 'cbind' method for \"DNAbin\" accepts only matrices")
115 NR <- unlist(lapply(obj, nrow))
116 for (i in 1:n) class(obj[[i]]) <- NULL
118 nms <- unlist(lapply(obj, rownames))
119 if (fill.with.gaps) {
120 NC <- unlist(lapply(obj, ncol))
122 ans <- matrix(as.raw(4), length(nms), sum(NC))
126 to <- from + NC[i] - 1
127 tmp <- rownames(obj[[i]])
128 nmsi <- tmp[tmp %in% nms]
129 ans[nmsi, from:to] <- obj[[i]][nmsi, , drop = FALSE]
135 nms <- names(tab)[which(ubi)]
136 ans <- obj[[1]][nms, , drop = FALSE]
138 ans <- cbind(ans, obj[[i]][nms, , drop = FALSE])
139 if (!quiet && !all(ubi))
140 warning("some rows were dropped.")
143 if (length(unique(NR)) > 1)
144 stop("matrices do not have the same number of rows.")
145 ans <- matrix(unlist(obj), NR)
146 rownames(ans) <- rownames(obj[[1]])
148 class(ans) <- "DNAbin"
152 c.DNAbin <- function(..., recursive = FALSE)
154 if (!all(unlist(lapply(list(...), is.list))))
155 stop("the 'c' method for \"DNAbin\" accepts only lists")
156 structure(NextMethod("c"), class = "DNAbin")
159 print.DNAbin <- function(x, ...)
161 n <- 1 # <- if is.vector(x)
162 if (is.list(x)) n <- length(x)
163 else if (is.matrix(x)) n <- dim(x)[1]
164 if (n > 1) cat(n, "DNA sequences in binary format.\n")
165 else cat("1 DNA sequence in binary format.\n")
168 summary.DNAbin <- function(object, printlen = 6, digits = 3, ...)
170 if (is.list(object)) {
174 cat("1 DNA sequence in binary format stored in a list.\n\n")
175 cat("Sequence length:", length(object[[1]]), "\n\n")
176 cat("Label:", nms, "\n\n")
178 cat(n, "DNA sequences in binary format stored in a list.\n\n")
179 tmp <- unlist(lapply(object, length))
183 cat("All sequences of same length:", maxi, "\n")
185 cat("Mean sequence length:", round(mean(tmp), 3), "\n")
186 cat(" Shortest sequence:", mini, "\n")
187 cat(" Longest sequence:", maxi, "\n")
191 nms <- nms[1:printlen]
194 cat("\nLabels:", paste(nms, collapse = " "), TAIL)
196 } else if (is.matrix(object)) {
198 nms <- rownames(object)
199 cat(nd[1], "DNA sequences in binary format stored in a matrix.\n\n")
200 cat("All sequences of same length:", nd[2], "\n")
202 if (printlen < nd[1]) {
203 nms <- nms[1:printlen]
206 cat("\nLabels:", paste(nms, collapse = " "), TAIL)
208 cat("1 DNA sequence in binary format stored in a vector.\n\n")
209 cat("Sequence length:", length(object), "\n\n")
211 cat("Base composition:\n")
212 print(round(base.freq(object), digits))
215 as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
217 ._cs_<- letters[c(1, 7, 3, 20, 18, 13, 23, 19, 11, 25, 22, 8, 4, 2, 14)]
219 ._bs_<- c(136, 72, 40, 24, 192, 160, 144, 96, 80, 48, 224, 176, 208, 112, 240)
221 as.DNAbin.character <- function(x, ...)
226 ans[which(x == ._cs_[i])] <- as.raw(._bs_[i])
227 ans[which(x == "-")] <- as.raw(4)
228 ans[which(x == "?")] <- as.raw(2)
231 dimnames(ans) <- dimnames(x)
233 class(ans) <- "DNAbin"
237 as.DNAbin.alignment <- function(x, ...)
240 x$seq <- tolower(x$seq)
241 ans <- matrix("", n, nchar(x$seq[1]))
243 ans[i, ] <- strsplit(x$seq[i], "")[[1]]
244 rownames(ans) <- gsub(" +$", "", gsub("^ +", "", x$nam))
245 as.DNAbin.character(ans)
248 as.DNAbin.list <- function(x, ...)
250 obj <- lapply(x, as.DNAbin)
251 class(obj) <- "DNAbin"
255 as.character.DNAbin <- function(x, ...)
258 ans <- character(length(xx))
260 ans[which(xx == ._bs_[i])] <- ._cs_[i]
261 ans[which(xx == 4)] <- "-"
262 ans[which(xx == 2)] <- "?"
265 dimnames(ans) <- dimnames(xx)
269 if (is.list(x)) lapply(x, f) else f(x)
272 base.freq <- function(x, freq = FALSE)
274 if (is.list(x)) x <- unlist(x)
276 BF <- .C("BaseProportion", x, n, double(4), freq,
277 DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
278 names(BF) <- letters[c(1, 3, 7, 20)]
282 GC.content <- function(x) sum(base.freq(x)[2:3])
284 seg.sites <- function(x)
286 if (is.list(x)) x <- as.matrix(x)
287 if (is.vector(x)) n <- 1
288 else { # 'x' is a matrix
293 if (n == 1) return(integer(0))
294 ans <- .C("SegSites", x, n, s, integer(s),
295 DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
296 which(as.logical(ans[[4]]))
299 dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
300 pairwise.deletion = FALSE, base.freq = NULL,
303 MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93",
304 "GG95", "LOGDET", "BH87", "PARALIN", "N")
305 imod <- pmatch(toupper(model), MODELS)
307 stop(paste("'model' must be one of:",
308 paste("\"", MODELS, "\"", sep = "", collapse = " ")))
309 if (imod == 11 && variance) {
310 warning("computing variance temporarily not available for model BH87.")
313 if (gamma && imod %in% c(1, 5:7, 9:12)) {
314 warning(paste("gamma-correction not available for model", model))
317 if (is.list(x)) x <- as.matrix(x)
318 nms <- dimnames(x)[[1]]
322 BF <- if (is.null(base.freq)) base.freq(x) else base.freq
323 if (!pairwise.deletion) {
324 keep <- .C("GlobalDeletionDNA", x, n, s,
325 rep(1L, s), PACKAGE = "ape")[[4]]
326 x <- x[, as.logical(keep)]
329 Ndist <- if (imod == 11) n*n else n*(n - 1)/2
330 var <- if (variance) double(Ndist) else 0
331 if (!gamma) gamma <- alpha <- 0
332 else alpha <- gamma <- 1
333 d <- .C("dist_dna", x, n, s, imod, double(Ndist), BF,
334 as.integer(pairwise.deletion), as.integer(variance),
335 var, as.integer(gamma), alpha, DUP = FALSE, NAOK = TRUE,
337 if (variance) var <- d[[9]]
341 dimnames(d) <- list(nms, nms)
344 attr(d, "Labels") <- nms
345 attr(d, "Diag") <- attr(d, "Upper") <- FALSE
346 attr(d, "call") <- match.call()
347 attr(d, "method") <- model
349 if (as.matrix) d <- as.matrix(d)
351 if (variance) attr(d, "variance") <- var