R/DNA.R

   1 ## DNA.R (2011-03-15)
   2
   3 ##   Manipulations and Comparisons of DNA Sequences
   4
   5 ## Copyright 2002-2011 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 as.list.DNAbin <- function(x, ...)
  95 {
  96     if (is.list(x)) return(x)
  97     if (is.null(dim(x))) obj <- list(x) # cause is.vector() doesn't work
  98     else { # matrix
  99         n <- nrow(x)
 100         obj <- vector("list", n)
 101         for (i in 1:n) obj[[i]] <- x[i, ]
 102         names(obj) <- rownames(x)
 103     }
 104     class(obj) <- "DNAbin"
 105     obj
 106 }
 107
 108 rbind.DNAbin <- function(...)
 109 ### works only with matrices for the moment
 110 {
 111     obj <- list(...)
 112     n <- length(obj)
 113     if (n == 1) return(obj[[1]])
 114     for (i in 1:n)
 115         if (!is.matrix(obj[[1]]))
 116             stop("the 'rbind' method for \"DNAbin\" accepts only matrices")
 117     NC <- unlist(lapply(obj, ncol))
 118     if (length(unique(NC)) > 1)
 119         stop("matrices do not have the same number of columns.")
 120     for (i in 1:n) class(obj[[i]]) <- NULL
 121     for (i in 2:n) obj[[1]] <- rbind(obj[[1]], obj[[i]])
 122     structure(obj[[1]], class = "DNAbin")
 123 }
 124
 125 cbind.DNAbin <-
 126     function(..., check.names = TRUE, fill.with.gaps = FALSE,
 127              quiet = FALSE)
 128 ### works only with matrices for the moment
 129 {
 130     obj <- list(...)
 131     n <- length(obj)
 132     if (n == 1) return(obj[[1]])
 133     for (i in 1:n)
 134         if (!is.matrix(obj[[1]]))
 135             stop("the 'cbind' method for \"DNAbin\" accepts only matrices")
 136     NR <- unlist(lapply(obj, nrow))
 137     for (i in 1:n) class(obj[[i]]) <- NULL
 138     if (check.names) {
 139         nms <- unlist(lapply(obj, rownames))
 140         if (fill.with.gaps) {
 141             NC <- unlist(lapply(obj, ncol))
 142             nms <- unique(nms)
 143             ans <- matrix(as.raw(4), length(nms), sum(NC))
 144             rownames(ans) <- nms
 145             from <- 1
 146             for (i in 1:n) {
 147                 to <- from + NC[i] - 1
 148                 tmp <- rownames(obj[[i]])
 149                 nmsi <- tmp[tmp %in% nms]
 150                 ans[nmsi, from:to] <- obj[[i]][nmsi, , drop = FALSE]
 151                 from <- to + 1
 152             }
 153         } else {
 154             tab <- table(nms)
 155             ubi <- tab == n
 156             nms <- names(tab)[which(ubi)]
 157             ans <- obj[[1]][nms, , drop = FALSE]
 158             for (i in 2:n)
 159                 ans <- cbind(ans, obj[[i]][nms, , drop = FALSE])
 160             if (!quiet && !all(ubi))
 161                 warning("some rows were dropped.")
 162         }
 163     } else {
 164         if (length(unique(NR)) > 1)
 165             stop("matrices do not have the same number of rows.")
 166         ans <- matrix(unlist(obj), NR)
 167         rownames(ans) <- rownames(obj[[1]])
 168     }
 169     class(ans) <- "DNAbin"
 170     ans
 171 }
 172
 173 c.DNAbin <- function(..., recursive = FALSE)
 174 {
 175     if (!all(unlist(lapply(list(...), is.list))))
 176         stop("the 'c' method for \"DNAbin\" accepts only lists")
 177     structure(NextMethod("c"), class = "DNAbin")
 178 }
 179
 180 print.DNAbin <- function(x, printlen = 6, digits = 3, ...)
 181 {
 182     if (is.list(x)) {
 183         n <- length(x)
 184         nms <- names(x)
 185         if (n == 1) {
 186             cat("1 DNA sequence in binary format stored in a list.\n\n")
 187             cat("Sequence length:", length(x[[1]]), "\n\n")
 188             cat("Label:", nms, "\n\n")
 189         } else {
 190             cat(n, "DNA sequences in binary format stored in a list.\n\n")
 191             tmp <- unlist(lapply(x, length))
 192             mini <- min(tmp)
 193             maxi <- max(tmp)
 194             if (mini == maxi)
 195                 cat("All sequences of same length:", maxi, "\n")
 196             else {
 197                 cat("Mean sequence length:", round(mean(tmp), 3), "\n")
 198                 cat("   Shortest sequence:", mini, "\n")
 199                 cat("    Longest sequence:", maxi, "\n")
 200             }
 201             TAIL <- "\n\n"
 202             if (printlen < n) {
 203                 nms <- nms[1:printlen]
 204                 TAIL <- "...\n\n"
 205             }
 206             cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 207         }
 208     } else if (is.matrix(x)) {
 209         nd <- dim(x)
 210         nms <- rownames(x)
 211         cat(nd[1], "DNA sequences in binary format stored in a matrix.\n\n")
 212         cat("All sequences of same length:", nd[2], "\n")
 213         TAIL <- "\n\n"
 214         if (printlen < nd[1]) {
 215             nms <- nms[1:printlen]
 216             TAIL <- "...\n\n"
 217         }
 218         cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 219     } else {
 220         cat("1 DNA sequence in binary format stored in a vector.\n\n")
 221         cat("Sequence length:", length(x), "\n\n")
 222     }
 223     cat("Base composition:\n")
 224     print(round(base.freq(x), digits))
 225 }
 226
 227 as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
 228
 229 ._cs_ <- c("a", "g", "c", "t", "r", "m", "w", "s", "k",
 230            "y", "v", "h",  "d", "b", "n", "-", "?")
 231
 232 ._bs_ <- c(136, 72, 40, 24, 192, 160, 144, 96, 80,
 233            48, 224, 176, 208, 112, 240, 4, 2)
 234
 235 as.DNAbin.character <- function(x, ...)
 236 {
 237     n <- length(x)
 238     ans <- raw(n)
 239     for (i in 1:15)
 240       ans[which(x == ._cs_[i])] <- as.raw(._bs_[i])
 241     ans[which(x == "-")] <- as.raw(4)
 242     ans[which(x == "?")] <- as.raw(2)
 243     if (is.matrix(x)) {
 244         dim(ans) <- dim(x)
 245         dimnames(ans) <- dimnames(x)
 246     }
 247     class(ans) <- "DNAbin"
 248     ans
 249 }
 250
 251 as.DNAbin.alignment <- function(x, ...)
 252 {
 253     n <- x$nb
 254     x$seq <- tolower(x$seq)
 255     ans <- matrix("", n, nchar(x$seq[1]))
 256     for (i in 1:n)
 257         ans[i, ] <- strsplit(x$seq[i], "")[[1]]
 258     rownames(ans) <- gsub(" +$", "", gsub("^ +", "", x$nam))
 259     as.DNAbin.character(ans)
 260 }
 261
 262 as.DNAbin.list <- function(x, ...)
 263 {
 264     obj <- lapply(x, as.DNAbin)
 265     class(obj) <- "DNAbin"
 266     obj
 267 }
 268
 269 as.character.DNAbin <- function(x, ...)
 270 {
 271     f <- function(xx) {
 272         ans <- character(length(xx))
 273         for (i in 1:15)
 274           ans[which(xx == ._bs_[i])] <- ._cs_[i]
 275         ans[which(xx == 4)] <- "-"
 276         ans[which(xx == 2)] <- "?"
 277         if (is.matrix(xx)) {
 278             dim(ans) <- dim(xx)
 279             dimnames(ans) <- dimnames(xx)
 280         }
 281         ans
 282     }
 283     if (is.list(x)) lapply(x, f) else f(x)
 284 }
 285
 286 base.freq <- function(x, freq = FALSE, all = FALSE)
 287 {
 288     if (is.list(x)) x <- unlist(x)
 289     n <- length(x)
 290     BF <-.C("BaseProportion", x, n, double(17),
 291             DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]]
 292     names(BF) <- c("a", "c", "g", "t", "r", "m", "w", "s",
 293                    "k", "y", "v", "h", "d", "b", "n", "-", "?")
 294     if (all) {
 295         if (!freq) BF <- BF / n
 296     } else {
 297         BF <- BF[1:4]
 298         if (!freq) BF <- BF / sum(BF)
 299     }
 300     BF
 301 }
 302
 303 Ftab <- function(x, y = NULL)
 304 {
 305     if (is.null(y)) {
 306         if (is.list(x)) {
 307             y <- x[[2]]
 308             x <- x[[1]]
 309             if (length(x) != length(y))
 310                 stop("'x' and 'y' not of same lenght")
 311         } else { # 'x' is a matrix
 312             y <- x[2, , drop = TRUE]
 313             x <- x[1, , drop = TRUE]
 314         }
 315     } else {
 316         x <- as.vector(x)
 317         y <- as.vector(y)
 318         if (length(x) != length(y))
 319             stop("'x' and 'y' not of same lenght")
 320     }
 321     out <- matrix(0, 4, 4)
 322     k <- c(136, 40, 72, 24)
 323     for (i in 1:4) {
 324         a <- x == k[i]
 325         for (j in 1:4) {
 326             b <- y == k[j]
 327             out[i, j] <- sum(a & b)
 328         }
 329     }
 330     dimnames(out)[1:2] <- list(c("a", "c", "g", "t"))
 331     out
 332 }
 333
 334 GC.content <- function(x) sum(base.freq(x)[2:3])
 335
 336 seg.sites <- function(x)
 337 {
 338     if (is.list(x)) x <- as.matrix(x)
 339     if (is.vector(x)) n <- 1
 340     else { # 'x' is a matrix
 341         n <- dim(x)
 342         s <- n[2]
 343         n <- n[1]
 344     }
 345     if (n == 1) return(integer(0))
 346     ans <- .C("SegSites", x, n, s, integer(s),
 347               DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
 348     which(as.logical(ans[[4]]))
 349 }
 350
 351 dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
 352                      pairwise.deletion = FALSE, base.freq = NULL,
 353                      as.matrix = FALSE)
 354 {
 355     MODELS <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93",
 356                 "GG95", "LOGDET", "BH87", "PARALIN", "N", "TS", "TV")
 357     imod <- pmatch(toupper(model), MODELS)
 358     if (is.na(imod))
 359         stop(paste("'model' must be one of:",
 360                    paste("\"", MODELS, "\"", sep = "", collapse = " ")))
 361     if (imod == 11 && variance) {
 362         warning("computing variance temporarily not available for model BH87.")
 363         variance <- FALSE
 364     }
 365     if (gamma && imod %in% c(1, 5:7, 9:15)) {
 366         warning(paste("gamma-correction not available for model", model))
 367         gamma <- FALSE
 368     }
 369     if (is.list(x)) x <- as.matrix(x)
 370     nms <- dimnames(x)[[1]]
 371     n <- dim(x)
 372     s <- n[2]
 373     n <- n[1]
 374     if (imod %in% c(4, 6:8)) {
 375         BF <- if (is.null(base.freq)) base.freq(x) else base.freq
 376     } else BF <- 0
 377     if (!pairwise.deletion) {
 378         keep <- .C("GlobalDeletionDNA", x, n, s,
 379                    rep(1L, s), PACKAGE = "ape")[[4]]
 380         x <- x[,  as.logical(keep)]
 381         s <- dim(x)[2]
 382     }
 383     Ndist <- if (imod == 11) n*n else n*(n - 1)/2
 384     var <- if (variance) double(Ndist) else 0
 385     if (!gamma) gamma <- alpha <- 0
 386     else alpha <- gamma <- 1
 387     d <- .C("dist_dna", x, n, s, imod, double(Ndist), BF,
 388             as.integer(pairwise.deletion), as.integer(variance),
 389             var, as.integer(gamma), alpha, DUP = FALSE, NAOK = TRUE,
 390             PACKAGE = "ape")
 391     if (variance) var <- d[[9]]
 392     d <- d[[5]]
 393     if (imod == 11) {
 394         dim(d) <- c(n, n)
 395         dimnames(d) <- list(nms, nms)
 396     } else {
 397         attr(d, "Size") <- n
 398         attr(d, "Labels") <- nms
 399         attr(d, "Diag") <- attr(d, "Upper") <- FALSE
 400         attr(d, "call") <- match.call()
 401         attr(d, "method") <- model
 402         class(d) <- "dist"
 403         if (as.matrix) d <- as.matrix(d)
 404     }
 405     if (variance) attr(d, "variance") <- var
 406     d
 407 }
 408
 409 image.DNAbin <- function(x, what, col, bg = "white", xlab = "", ylab = "",
 410                          show.labels = TRUE, cex.lab = 1, legend = TRUE, ...)
 411 {
 412     what <-
 413         if (missing(what)) c("a", "g", "c", "t", "n", "-") else tolower(what)
 414     if (missing(col))
 415         col <- c("red", "yellow", "green", "blue", "grey", "black")
 416     n <- (dx <- dim(x))[1] # number of sequences
 417     s <- dx[2] # number of sites
 418     y <- integer(N <- length(x))
 419     ncl <- length(what)
 420     col <- rep(col, length.out = ncl)
 421     sm <- 0L
 422     for (i in ncl:1) {
 423         k <- ._bs_[._cs_ == what[i]]
 424         sel <- which(x == k)
 425         if (ll <- length(sel)) {
 426             y[sel] <- i
 427             sm <- sm + ll
 428         } else {
 429             what <- what[-i]
 430             col <- col[-i]
 431         }
 432     }
 433     dim(y) <- dx
 434     ## if there's no 0 in y, must drop 'bg' from the cols passed to image:
 435     if (sm == N) {
 436         leg.co <- co <- col
 437         leg.txt <- toupper(what)
 438     } else {
 439         co <- c(bg, col)
 440         leg.txt <- c(toupper(what), "others")
 441         leg.co <- c(col, bg)
 442     }
 443     yaxt <- if (show.labels) "n" else "s"
 444     image(1:s, 1:n, t(y), col = co, xlab = xlab,
 445           ylab = ylab, yaxt = yaxt, ...)
 446     if (show.labels)
 447         mtext(rownames(x), side = 2, line = 0.1, at = 1:n,
 448               cex = cex.lab, adj = 1, las = 1)
 449     if (legend) {
 450         psr <- par("usr")
 451         xx <- psr[2]/2
 452         yy <- psr[4] * (0.5 + 0.5/par("plt")[4])
 453         legend(xx, yy, legend = leg.txt, pch = 22, pt.bg = leg.co,
 454                pt.cex = 2, bty = "n", xjust = 0.5, yjust = 0.5,
 455                horiz = TRUE, xpd = TRUE)
 456     }
 457 }