R/DNA.R

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