R/DNA.R

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