3 ## Manipulations and Comparisons of DNA Sequences
5 ## Copyright 2002-2008 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 as.alignment <- function(x)
12 if (is.list(x)) n <- length(x)
13 if (is.matrix(x)) n <- dim(x)[1]
18 seq[i] <- paste(x[[i]], collapse = "")
21 nam <- dimnames(x)[[1]]
23 seq[i] <- paste(x[i, ], collapse = "")
25 obj <- list(nb = n, seq = seq, nam = nam, com = NA)
26 class(obj) <- "alignment"
30 "[.DNAbin" <- function(x, i, j, drop = TRUE)
34 if (nargs() == 2 && !missing(i)) ans <- x[i]
37 if (missing(i)) i <- 1:nd[1]
38 if (missing(j)) j <- 1:nd[2]
39 ans <- x[i, j, drop = drop]
42 if (missing(i)) i <- 1:length(x)
45 structure(ans, class = "DNAbin")
48 as.matrix.DNAbin <- function(x, ...)
50 if (is.matrix(x)) return(x)
52 if (length(unique(unlist(lapply(x, length)))) != 1)
53 stop("DNA sequences in list not of the same length.")
57 x <- matrix(unlist(x), n, s, byrow = TRUE)
64 rbind.DNAbin <- function(...)
65 ### works only with matrices for the moment
69 if (nobj == 1) stop("only one matrix to bind.")
72 if(ncol(obj[[i]]) != NC)
73 stop("matrices do not have the same number of columns.")
74 for (i in 1:nobj) class(obj[[i]]) <- NULL
76 for (i in 2:nobj) ans <- rbind(ans, obj[[i]])
77 structure(ans, class = "DNAbin")
80 cbind.DNAbin <- function(..., check.names = TRUE)
81 ### works only with matrices for the moment
85 if (nobj == 1) stop("only one matrix to bind.")
88 if(nrow(obj[[i]]) != NR)
89 stop("matrices do not have the same number of rows.")
90 for (i in 1:nobj) class(obj[[i]]) <- NULL
91 nms <- rownames(obj[[1]])
94 if (all(rownames(obj[[i]]) %in% nms))
95 obj[[i]] <- obj[[i]][nms, ]
96 else stop("rownames do not match among matrices.")
98 ans <- matrix(unlist(obj), NR)
100 structure(ans, class = "DNAbin")
103 print.DNAbin <- function(x, ...)
105 n <- 1 # <- if is.vector(x)
106 if (is.list(x)) n <- length(x)
107 else if (is.matrix(x)) n <- dim(x)[1]
108 if (n > 1) cat(n, "DNA sequences in binary format.\n")
109 else cat("1 DNA sequence in binary format.\n")
112 summary.DNAbin <- function(object, printlen = 6, digits = 3, ...)
114 if (is.list(object)) {
118 cat("1 DNA sequence in binary format stored in a list.\n\n")
119 cat("Sequence length:", length(object[[1]]), "\n\n")
120 cat("Label:", nms, "\n\n")
122 cat(n, "DNA sequences in binary format stored in a list.\n\n")
123 cat("Summary of sequence lengths:\n")
124 print(summary(unlist(lapply(object, length))))
127 nms <- nms[1:printlen]
130 cat("\nLabels:", paste(nms, collapse = " "), TAIL)
132 } else if (is.matrix(object)) {
134 nms <- rownames(object)
135 cat(nd[1], "DNA sequences in binary format stored in a matrix.\n\n")
136 cat("All sequences of same length:", nd[2], "\n")
138 if (printlen < nd[1]) {
139 nms <- nms[1:printlen]
142 cat("\nLabels:", paste(nms, collapse = " "), TAIL)
144 cat("1 DNA sequence in binary format stored in a vector.\n\n")
145 cat("Sequence length:", length(object), "\n\n")
147 cat("Base composition:\n")
148 print(round(base.freq(object), digits))
151 as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
153 ._cs_<- letters[c(1, 7, 3, 20, 18, 13, 23, 19, 11, 25, 22, 8, 4, 2, 14)]
155 ._bs_<- c(136, 72, 40, 24, 192, 160, 144, 96, 80, 48, 224, 176, 208, 112, 240)
157 as.DNAbin.character <- function(x, ...)
162 ans[which(x == ._cs_[i])] <- as.raw(._bs_[i])
163 ans[which(x == "-")] <- as.raw(4)
164 ans[which(x == "?")] <- as.raw(2)
167 dimnames(ans) <- dimnames(x)
169 class(ans) <- "DNAbin"
173 as.DNAbin.alignment <- function(x, ...)
176 x$seq <- tolower(x$seq)
177 ans <- matrix("", n, nchar(x$seq[1]))
179 ans[i, ] <- strsplit(x$seq[i], "")[[1]]
180 rownames(ans) <- gsub(" +$", "", gsub("^ +", "", x$nam))
181 as.DNAbin.character(ans)
184 as.DNAbin.list <- function(x, ...)
186 obj <- lapply(x, as.DNAbin)
187 class(obj) <- "DNAbin"
191 as.character.DNAbin <- function(x, ...)
194 ans <- character(length(xx))
196 ans[which(xx == ._bs_[i])] <- ._cs_[i]
197 ans[which(xx == 4)] <- "-"
198 ans[which(xx == 2)] <- "?"
201 dimnames(ans) <- dimnames(xx)
205 if (is.list(x)) lapply(x, f) else f(x)
208 base.freq <- function(x)
210 if (is.list(x)) x <- unlist(x)
212 BF <- .C("BaseProportion", as.raw(x), as.integer(n),
213 double(4), PACKAGE = "ape")[[3]]
214 names(BF) <- letters[c(1, 3, 7, 20)]
218 GC.content <- function(x)
224 seg.sites <- function(x)
229 ans <- .C("SegSites", x, as.integer(n), as.integer(s),
230 integer(s), PACKAGE = "ape")
231 which(as.logical(ans[[4]]))
234 nuc.div <- function(x, variance = FALSE, pairwise.deletion = FALSE)
236 if (pairwise.deletion && variance)
237 warning("cannot compute the variance of nucleotidic diversity\nwith pairwise deletion: try 'pairwise.deletion = FALSE' instead.")
238 if (is.list(x)) x <- as.matrix(x)
240 ans <- sum(dist.dna(x, "raw", pairwise.deletion = pairwise.deletion))/
243 var <- (n + 1)*ans/(3*(n + 1)*dim(x)[2]) + 2*(n^2 + n + 3)*ans/(9*n*(n - 1))
249 dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
250 pairwise.deletion = FALSE, base.freq = NULL,
253 MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93",
254 "GG95", "LOGDET", "BH87", "PARALIN")
255 imod <- which(MODELS == toupper(model))
256 if (imod == 11 && variance) {
257 warning("computing variance temporarily not available for model BH87.")
260 if (gamma && imod %in% c(1, 5:7, 9:12)) {
261 warning(paste("gamma-correction not available for model", model))
264 if (is.list(x)) x <- as.matrix(x)
265 nms <- dimnames(x)[[1]]
269 BF <- if (is.null(base.freq)) base.freq(x) else base.freq
270 if (!pairwise.deletion) {
271 keep <- .C("GlobalDeletionDNA", x, as.integer(n),
272 as.integer(s), as.integer(rep(1, s)),
273 PACKAGE = "ape")[[4]]
274 x <- x[, as.logical(keep)]
277 Ndist <- if (imod == 11) n*n else n*(n - 1)/2
278 var <- if (variance) double(Ndist) else 0
279 if (!gamma) gamma <- alpha <- 0
280 else alpha <- gamma <- 1
281 d <- .C("dist_dna", x, as.integer(n), as.integer(s),
282 as.integer(imod), double(Ndist), BF,
283 as.integer(pairwise.deletion), as.integer(variance),
284 var, as.integer(gamma), alpha, PACKAGE = "ape")
285 if (variance) var <- d[[9]]
289 dimnames(d) <- list(nms, nms)
292 attr(d, "Labels") <- nms
293 attr(d, "Diag") <- attr(d, "Upper") <- FALSE
294 attr(d, "call") <- match.call()
295 attr(d, "method") <- model
297 if (as.matrix) d <- as.matrix(d)
299 if (variance) attr(d, "variance") <- var