R/DNA.R

   1 ## DNA.R (2009-05-19)
   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 c.DNAbin <- function(..., recursive = FALSE)
 146     structure(NextMethod("c"), class = "DNAbin")
 147
 148 print.DNAbin <- function(x, ...)
 149 {
 150     n <- 1 # <- if is.vector(x)
 151     if (is.list(x)) n <- length(x)
 152     else if (is.matrix(x)) n <- dim(x)[1]
 153     if (n > 1) cat(n, "DNA sequences in binary format.\n")
 154     else cat("1 DNA sequence in binary format.\n")
 155 }
 156
 157 summary.DNAbin <- function(object, printlen = 6, digits = 3, ...)
 158 {
 159     if (is.list(object)) {
 160         n <- length(object)
 161         nms <- names(object)
 162         if (n == 1) {
 163             cat("1 DNA sequence in binary format stored in a list.\n\n")
 164             cat("Sequence length:", length(object[[1]]), "\n\n")
 165             cat("Label:", nms, "\n\n")
 166         } else {
 167             cat(n, "DNA sequences in binary format stored in a list.\n\n")
 168             tmp <- unlist(lapply(object, length))
 169             mini <- min(tmp)
 170             maxi <- max(tmp)
 171             if (mini == maxi)
 172                 cat("All sequences of same length:", maxi, "\n")
 173             else {
 174                 cat("Mean sequence length:", round(mean(tmp), 3), "\n")
 175                 cat("   Shortest sequence:", mini, "\n")
 176                 cat("    Longest sequence:", maxi, "\n")
 177             }
 178             TAIL <- "\n\n"
 179             if (printlen < n) {
 180                 nms <- nms[1:printlen]
 181                 TAIL <- "...\n\n"
 182             }
 183             cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 184         }
 185     } else if (is.matrix(object)) {
 186         nd <- dim(object)
 187         nms <- rownames(object)
 188         cat(nd[1], "DNA sequences in binary format stored in a matrix.\n\n")
 189         cat("All sequences of same length:", nd[2], "\n")
 190         TAIL <- "\n\n"
 191         if (printlen < nd[1]) {
 192             nms <- nms[1:printlen]
 193             TAIL <- "...\n\n"
 194         }
 195         cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 196     } else {
 197         cat("1 DNA sequence in binary format stored in a vector.\n\n")
 198         cat("Sequence length:", length(object), "\n\n")
 199     }
 200     cat("Base composition:\n")
 201     print(round(base.freq(object), digits))
 202 }
 203
 204 as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
 205
 206 ._cs_<- letters[c(1, 7, 3, 20, 18, 13, 23, 19, 11, 25, 22, 8, 4, 2, 14)]
 207
 208 ._bs_<- c(136, 72, 40, 24, 192, 160, 144, 96, 80, 48, 224, 176, 208, 112, 240)
 209
 210 as.DNAbin.character <- function(x, ...)
 211 {
 212     n <- length(x)
 213     ans <- raw(n)
 214     for (i in 1:15)
 215       ans[which(x == ._cs_[i])] <- as.raw(._bs_[i])
 216     ans[which(x == "-")] <- as.raw(4)
 217     ans[which(x == "?")] <- as.raw(2)
 218     if (is.matrix(x)) {
 219         dim(ans) <- dim(x)
 220         dimnames(ans) <- dimnames(x)
 221     }
 222     class(ans) <- "DNAbin"
 223     ans
 224 }
 225
 226 as.DNAbin.alignment <- function(x, ...)
 227 {
 228     n <- x$nb
 229     x$seq <- tolower(x$seq)
 230     ans <- matrix("", n, nchar(x$seq[1]))
 231     for (i in 1:n)
 232         ans[i, ] <- strsplit(x$seq[i], "")[[1]]
 233     rownames(ans) <- gsub(" +$", "", gsub("^ +", "", x$nam))
 234     as.DNAbin.character(ans)
 235 }
 236
 237 as.DNAbin.list <- function(x, ...)
 238 {
 239     obj <- lapply(x, as.DNAbin)
 240     class(obj) <- "DNAbin"
 241     obj
 242 }
 243
 244 as.character.DNAbin <- function(x, ...)
 245 {
 246     f <- function(xx) {
 247         ans <- character(length(xx))
 248         for (i in 1:15)
 249           ans[which(xx == ._bs_[i])] <- ._cs_[i]
 250         ans[which(xx == 4)] <- "-"
 251         ans[which(xx == 2)] <- "?"
 252         if (is.matrix(xx)) {
 253             dim(ans) <- dim(xx)
 254             dimnames(ans) <- dimnames(xx)
 255         }
 256         ans
 257     }
 258     if (is.list(x)) lapply(x, f) else f(x)
 259 }
 260
 261 base.freq <- function(x)
 262 {
 263     if (is.list(x)) x <- unlist(x)
 264     n <- length(x)
 265     BF <- .C("BaseProportion", x, n, double(4),
 266              DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
 267     names(BF) <- letters[c(1, 3, 7, 20)]
 268     BF
 269 }
 270
 271 GC.content <- function(x) sum(base.freq(x)[2:3])
 272
 273 seg.sites <- function(x)
 274 {
 275     if (is.list(x)) x <- as.matrix(x)
 276     n <- dim(x)
 277     s <- n[2]
 278     n <- n[1]
 279     ans <- .C("SegSites", x, n, s, integer(s),
 280               DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
 281     which(as.logical(ans[[4]]))
 282 }
 283
 284 nuc.div <- function(x, variance = FALSE, pairwise.deletion = FALSE)
 285 {
 286     if (pairwise.deletion && variance)
 287       warning("cannot compute the variance of nucleotidic diversity\nwith pairwise deletion: try 'pairwise.deletion = FALSE' instead.")
 288     if (is.list(x)) x <- as.matrix(x)
 289     n <- dim(x)[1]
 290     ans <- sum(dist.dna(x, "raw", pairwise.deletion = pairwise.deletion))/
 291         (n*(n - 1)/2)
 292     if (variance) {
 293         var <- (n + 1)*ans/(3*(n + 1)*dim(x)[2]) + 2*(n^2 + n + 3)*ans/(9*n*(n - 1))
 294         ans <- c(ans, var)
 295     }
 296     ans
 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 <- which(MODELS == toupper(model))
 306     if (imod == 11 && variance) {
 307         warning("computing variance temporarily not available for model BH87.")
 308         variance <- FALSE
 309     }
 310     if (gamma && imod %in% c(1, 5:7, 9:12)) {
 311         warning(paste("gamma-correction not available for model", model))
 312         gamma <- FALSE
 313     }
 314     if (is.list(x)) x <- as.matrix(x)
 315     nms <- dimnames(x)[[1]]
 316     n <- dim(x)
 317     s <- n[2]
 318     n <- n[1]
 319     BF <- if (is.null(base.freq)) base.freq(x) else base.freq
 320     if (!pairwise.deletion) {
 321         keep <- .C("GlobalDeletionDNA", x, n, s,
 322                    rep(1L, s), PACKAGE = "ape")[[4]]
 323         x <- x[,  as.logical(keep)]
 324         s <- dim(x)[2]
 325     }
 326     Ndist <- if (imod == 11) n*n else n*(n - 1)/2
 327     var <- if (variance) double(Ndist) else 0
 328     if (!gamma) gamma <- alpha <- 0
 329     else alpha <- gamma <- 1
 330     d <- .C("dist_dna", x, n, s, imod, double(Ndist), BF,
 331             as.integer(pairwise.deletion), as.integer(variance),
 332             var, as.integer(gamma), alpha, DUP = FALSE, NAOK = TRUE,
 333             PACKAGE = "ape")
 334     if (variance) var <- d[[9]]
 335     d <- d[[5]]
 336     if (imod == 11) {
 337         dim(d) <- c(n, n)
 338         dimnames(d) <- list(nms, nms)
 339     } else {
 340         attr(d, "Size") <- n
 341         attr(d, "Labels") <- nms
 342         attr(d, "Diag") <- attr(d, "Upper") <- FALSE
 343         attr(d, "call") <- match.call()
 344         attr(d, "method") <- model
 345         class(d) <- "dist"
 346         if (as.matrix) d <- as.matrix(d)
 347     }
 348     if (variance) attr(d, "variance") <- var
 349     d
 350 }