R/DNA.R

   1 ## DNA.R (2010-03-16)
   2
   3 ##   Manipulations and Comparisons of DNA Sequences
   4
   5 ## Copyright 2002-2010 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 del.gaps <- function(x)
  11 {
  12     deleteGaps <- function(x) {
  13         i <- which(x == 4)
  14         if (length(i)) x[-i] else x
  15     }
  16
  17     if (!inherits(x, "DNAbin")) x <- as.DNAbin(x)
  18     if (is.matrix(x)) {
  19         n <- dim(x)[1]
  20         y <- vector("list", n)
  21         for (i in 1:n) y[[i]] <- x[i, ]
  22         names(y) <- rownames(x)
  23         x <- y
  24         rm(y)
  25     }
  26     if (!is.list(x)) return(deleteGaps(x))
  27     x <- lapply(x, deleteGaps)
  28     class(x) <- "DNAbin"
  29     x
  30 }
  31
  32 as.alignment <- function(x)
  33 {
  34     if (is.list(x)) n <- length(x)
  35     if (is.matrix(x)) n <- dim(x)[1]
  36     seq <- character(n)
  37     if (is.list(x)) {
  38         nam <- names(x)
  39         for (i in 1:n)
  40           seq[i] <- paste(x[[i]], collapse = "")
  41     }
  42     if (is.matrix(x)) {
  43         nam <- dimnames(x)[[1]]
  44         for (i in 1:n)
  45           seq[i] <- paste(x[i, ], collapse = "")
  46     }
  47     obj <- list(nb = n, seq = seq, nam = nam, com = NA)
  48     class(obj) <- "alignment"
  49     obj
  50 }
  51
  52 "[.DNAbin" <- function(x, i, j, drop = TRUE)
  53 {
  54     oc <- oldClass(x)
  55     class(x) <- NULL
  56     if (is.matrix(x)) {
  57         if (nargs() == 2 && !missing(i)) ans <- x[i]
  58         else {
  59             nd <- dim(x)
  60             if (missing(i)) i <- 1:nd[1]
  61             if (missing(j)) j <- 1:nd[2]
  62             ans <- x[i, j, drop = drop]
  63         }
  64     } else {
  65         if (missing(i)) i <- 1:length(x)
  66         ans <- x[i]
  67     }
  68     class(ans) <- oc
  69     ans
  70 }
  71
  72 as.matrix.DNAbin <- function(x, ...)
  73 {
  74     if (is.list(x)) {
  75         if (length(unique(unlist(lapply(x, length)))) != 1)
  76           stop("DNA sequences in list not of the same length.")
  77         nms <- names(x)
  78         n <- length(x)
  79         s <- length(x[[1]])
  80         x <- matrix(unlist(x), n, s, byrow = TRUE)
  81         rownames(x) <- nms
  82         class(x) <- "DNAbin"
  83     }
  84     x
  85 }
  86
  87 rbind.DNAbin <- function(...)
  88 ### works only with matrices for the moment
  89 {
  90     obj <- list(...)
  91     n <- length(obj)
  92     if (n == 1) return(obj[[1]])
  93     for (i in 1:n)
  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")
 102 }
 103
 104 cbind.DNAbin <-
 105     function(..., check.names = TRUE, fill.with.gaps = FALSE,
 106              quiet = FALSE)
 107 ### works only with matrices for the moment
 108 {
 109     obj <- list(...)
 110     n <- length(obj)
 111     if (n == 1) return(obj[[1]])
 112     for (i in 1:n)
 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
 117     if (check.names) {
 118         nms <- unlist(lapply(obj, rownames))
 119         if (fill.with.gaps) {
 120             NC <- unlist(lapply(obj, ncol))
 121             nms <- unique(nms)
 122             ans <- matrix(as.raw(4), length(nms), sum(NC))
 123             rownames(ans) <- nms
 124             from <- 1
 125             for (i in 1:n) {
 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]
 130                 from <- to + 1
 131             }
 132         } else {
 133             tab <- table(nms)
 134             ubi <- tab == n
 135             nms <- names(tab)[which(ubi)]
 136             ans <- obj[[1]][nms, , drop = FALSE]
 137             for (i in 2:n)
 138                 ans <- cbind(ans, obj[[i]][nms, , drop = FALSE])
 139             if (!quiet && !all(ubi))
 140                 warning("some rows were dropped.")
 141         }
 142     } else {
 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]])
 147     }
 148     class(ans) <- "DNAbin"
 149     ans
 150 }
 151
 152 c.DNAbin <- function(..., recursive = FALSE)
 153 {
 154     if (!all(unlist(lapply(list(...), is.list))))
 155         stop("the 'c' method for \"DNAbin\" accepts only lists")
 156     structure(NextMethod("c"), class = "DNAbin")
 157 }
 158
 159 print.DNAbin <- function(x, ...)
 160 {
 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")
 166 }
 167
 168 summary.DNAbin <- function(object, printlen = 6, digits = 3, ...)
 169 {
 170     if (is.list(object)) {
 171         n <- length(object)
 172         nms <- names(object)
 173         if (n == 1) {
 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")
 177         } else {
 178             cat(n, "DNA sequences in binary format stored in a list.\n\n")
 179             tmp <- unlist(lapply(object, length))
 180             mini <- min(tmp)
 181             maxi <- max(tmp)
 182             if (mini == maxi)
 183                 cat("All sequences of same length:", maxi, "\n")
 184             else {
 185                 cat("Mean sequence length:", round(mean(tmp), 3), "\n")
 186                 cat("   Shortest sequence:", mini, "\n")
 187                 cat("    Longest sequence:", maxi, "\n")
 188             }
 189             TAIL <- "\n\n"
 190             if (printlen < n) {
 191                 nms <- nms[1:printlen]
 192                 TAIL <- "...\n\n"
 193             }
 194             cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 195         }
 196     } else if (is.matrix(object)) {
 197         nd <- dim(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")
 201         TAIL <- "\n\n"
 202         if (printlen < nd[1]) {
 203             nms <- nms[1:printlen]
 204             TAIL <- "...\n\n"
 205         }
 206         cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 207     } else {
 208         cat("1 DNA sequence in binary format stored in a vector.\n\n")
 209         cat("Sequence length:", length(object), "\n\n")
 210     }
 211     cat("Base composition:\n")
 212     print(round(base.freq(object), digits))
 213 }
 214
 215 as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
 216
 217 ._cs_<- letters[c(1, 7, 3, 20, 18, 13, 23, 19, 11, 25, 22, 8, 4, 2, 14)]
 218
 219 ._bs_<- c(136, 72, 40, 24, 192, 160, 144, 96, 80, 48, 224, 176, 208, 112, 240)
 220
 221 as.DNAbin.character <- function(x, ...)
 222 {
 223     n <- length(x)
 224     ans <- raw(n)
 225     for (i in 1:15)
 226       ans[which(x == ._cs_[i])] <- as.raw(._bs_[i])
 227     ans[which(x == "-")] <- as.raw(4)
 228     ans[which(x == "?")] <- as.raw(2)
 229     if (is.matrix(x)) {
 230         dim(ans) <- dim(x)
 231         dimnames(ans) <- dimnames(x)
 232     }
 233     class(ans) <- "DNAbin"
 234     ans
 235 }
 236
 237 as.DNAbin.alignment <- function(x, ...)
 238 {
 239     n <- x$nb
 240     x$seq <- tolower(x$seq)
 241     ans <- matrix("", n, nchar(x$seq[1]))
 242     for (i in 1:n)
 243         ans[i, ] <- strsplit(x$seq[i], "")[[1]]
 244     rownames(ans) <- gsub(" +$", "", gsub("^ +", "", x$nam))
 245     as.DNAbin.character(ans)
 246 }
 247
 248 as.DNAbin.list <- function(x, ...)
 249 {
 250     obj <- lapply(x, as.DNAbin)
 251     class(obj) <- "DNAbin"
 252     obj
 253 }
 254
 255 as.character.DNAbin <- function(x, ...)
 256 {
 257     f <- function(xx) {
 258         ans <- character(length(xx))
 259         for (i in 1:15)
 260           ans[which(xx == ._bs_[i])] <- ._cs_[i]
 261         ans[which(xx == 4)] <- "-"
 262         ans[which(xx == 2)] <- "?"
 263         if (is.matrix(xx)) {
 264             dim(ans) <- dim(xx)
 265             dimnames(ans) <- dimnames(xx)
 266         }
 267         ans
 268     }
 269     if (is.list(x)) lapply(x, f) else f(x)
 270 }
 271
 272 base.freq <- function(x, freq = FALSE)
 273 {
 274     if (is.list(x)) x <- unlist(x)
 275     n <- length(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)]
 279     BF
 280 }
 281
 282 GC.content <- function(x) sum(base.freq(x)[2:3])
 283
 284 seg.sites <- function(x)
 285 {
 286     if (is.list(x)) x <- as.matrix(x)
 287     if (is.vector(x)) n <- 1
 288     else { # 'x' is a matrix
 289         n <- dim(x)
 290         s <- n[2]
 291         n <- n[1]
 292     }
 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]]))
 297 }
 298
 299 dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
 300                      pairwise.deletion = FALSE, base.freq = NULL,
 301                      as.matrix = FALSE)
 302 {
 303     MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93",
 304                 "GG95", "LOGDET", "BH87", "PARALIN", "N")
 305     imod <- pmatch(toupper(model), MODELS)
 306     if (is.na(imod))
 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.")
 311         variance <- FALSE
 312     }
 313     if (gamma && imod %in% c(1, 5:7, 9:12)) {
 314         warning(paste("gamma-correction not available for model", model))
 315         gamma <- FALSE
 316     }
 317     if (is.list(x)) x <- as.matrix(x)
 318     nms <- dimnames(x)[[1]]
 319     n <- dim(x)
 320     s <- n[2]
 321     n <- n[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)]
 327         s <- dim(x)[2]
 328     }
 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,
 336             PACKAGE = "ape")
 337     if (variance) var <- d[[9]]
 338     d <- d[[5]]
 339     if (imod == 11) {
 340         dim(d) <- c(n, n)
 341         dimnames(d) <- list(nms, nms)
 342     } else {
 343         attr(d, "Size") <- n
 344         attr(d, "Labels") <- nms
 345         attr(d, "Diag") <- attr(d, "Upper") <- FALSE
 346         attr(d, "call") <- match.call()
 347         attr(d, "method") <- model
 348         class(d) <- "dist"
 349         if (as.matrix) d <- as.matrix(d)
 350     }
 351     if (variance) attr(d, "variance") <- var
 352     d
 353 }