R/DNA.R

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