R/DNA.R

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