R/DNA.R

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