R/protein.R

   1 #' [.proteinbin
   2 #'
   3 #'
   4 "[.proteinbin" <- function(x,i,j,drop=FALSE) {
   5     ans <- NextMethod("[",drop=drop)
   6     class(ans) <- "proteinbin"
   7     ans
   8 }
   9
  10 as.matrix.proteinbin <- function(x,...) {
  11     if (is.matrix(x)) return(x)
  12     if (is.vector(x)) {
  13         dim(x) <- c(1, length(x))
  14         return(x)
  15     }
  16     s <- unique(unlist(lapply(x, length)))
  17     if (length(s) != 1)
  18         stop("Protein sequences in list not of the same length.")
  19     nms <- names(x)
  20     n <- length(x)
  21     y <- matrix(as.raw(0), n, s)
  22     for (i in seq_len(n)) y[i, ] <- x[[i]]
  23     rownames(y) <- nms
  24     class(y) <- "proteinbin"
  25     y
  26 }
  27
  28 as.list.proteinbin <- function(x, ...)
  29 {
  30     if (is.list(x)) return(x)
  31     if (is.null(dim(x))) obj <- list(x) # cause is.vector() doesn't work
  32     else { # matrix
  33         n <- nrow(x)
  34         obj <- vector("list", n)
  35         for (i in 1:n) obj[[i]] <- x[i, ]
  36         names(obj) <- rownames(x)
  37     }
  38     class(obj) <- "proteinbin"
  39     obj
  40 }
  41
  42 rbind.proteinbin <- function(...)
  43 ### works only with matrices for the moment
  44 {
  45     obj <- list(...)
  46     n <- length(obj)
  47     if (n == 1) return(obj[[1]])
  48     for (i in 1:n)
  49         if (!is.matrix(obj[[i]]))
  50             stop("the 'rbind' method for \"proteinbin\" accepts only matrices")
  51     NC <- unlist(lapply(obj, ncol))
  52     if (length(unique(NC)) > 1)
  53         stop("matrices do not have the same number of columns.")
  54     for (i in 1:n) class(obj[[i]]) <- NULL
  55     for (i in 2:n) obj[[1]] <- rbind(obj[[1]], obj[[i]])
  56     structure(obj[[1]], class = "proteinbin")
  57 }
  58
  59 cbind.proteinbin <-
  60     function(..., check.names = TRUE, fill.with.gaps = FALSE,
  61              quiet = FALSE)
  62 ### works only with matrices for the moment
  63 {
  64     obj <- list(...)
  65     n <- length(obj)
  66     if (n == 1) return(obj[[1]])
  67     for (i in 1:n)
  68         if (!is.matrix(obj[[i]]))
  69             stop("the 'cbind' method for \"proteinbin\" accepts only matrices")
  70     NR <- unlist(lapply(obj, nrow))
  71     for (i in 1:n) class(obj[[i]]) <- NULL
  72     if (check.names) {
  73         nms <- unlist(lapply(obj, rownames))
  74         if (fill.with.gaps) {
  75             NC <- unlist(lapply(obj, ncol))
  76             nms <- unique(nms)
  77             ans <- matrix(as.raw(4), length(nms), sum(NC))
  78             rownames(ans) <- nms
  79             from <- 1
  80             for (i in 1:n) {
  81                 to <- from + NC[i] - 1
  82                 tmp <- rownames(obj[[i]])
  83                 nmsi <- tmp[tmp %in% nms]
  84                 ans[nmsi, from:to] <- obj[[i]][nmsi, , drop = FALSE]
  85                 from <- to + 1
  86             }
  87         } else {
  88             tab <- table(nms)
  89             ubi <- tab == n
  90             nms <- names(tab)[which(ubi)]
  91             ans <- obj[[1]][nms, , drop = FALSE]
  92             for (i in 2:n)
  93                 ans <- cbind(ans, obj[[i]][nms, , drop = FALSE])
  94             if (!quiet && !all(ubi))
  95                 warning("some rows were dropped.")
  96         }
  97     } else {
  98         if (length(unique(NR)) > 1)
  99             stop("matrices do not have the same number of rows.")
 100         ans <- matrix(unlist(obj), NR)
 101         rownames(ans) <- rownames(obj[[1]])
 102     }
 103     class(ans) <- "proteinbin"
 104     ans
 105 }
 106
 107 c.proteinbin <- function(..., recursive = FALSE)
 108 {
 109     if (!all(unlist(lapply(list(...), is.list))))
 110         stop("the 'c' method for \"proteinbin\" accepts only lists")
 111     structure(NextMethod("c"), class = "proteinbin")
 112 }
 113
 114 print.proteinbin <- function(x, printlen = 6, digits = 3, ...)
 115 {
 116     if (is.list(x)) {
 117         n <- length(x)
 118         nms <- names(x)
 119         if (n == 1) {
 120             cat("1 Protein sequence stored in a list.\n\n")
 121             nTot <- length(x[[1]])
 122             cat("Sequence length:", nTot, "\n\n")
 123             cat("Label:", nms, "\n\n")
 124         } else {
 125             cat(n, "Protein sequences stored in a list.\n\n")
 126             tmp <- unlist(lapply(x, length))
 127             nTot <- sum(tmp)
 128             mini <- range(tmp)
 129             maxi <- mini[2]
 130             mini <- mini[1]
 131             if (mini == maxi)
 132                 cat("All sequences of same length:", maxi, "\n")
 133             else {
 134                 cat("Mean sequence length:", round(mean(tmp), 3), "\n")
 135                 cat("   Shortest sequence:", mini, "\n")
 136                 cat("    Longest sequence:", maxi, "\n")
 137             }
 138             TAIL <- "\n\n"
 139             if (printlen < n) {
 140                 nms <- nms[1:printlen]
 141                 TAIL <- "...\n\n"
 142             }
 143             cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 144         }
 145     } else {
 146         nTot <- length(x)
 147         if (is.matrix(x)) {
 148             nd <- dim(x)
 149             nms <- rownames(x)
 150             cat(nd[1], "Protein sequences stored in a matrix.\n\n")
 151             cat("All sequences of same length:", nd[2], "\n")
 152             TAIL <- "\n\n"
 153             if (printlen < nd[1]) {
 154                 nms <- nms[1:printlen]
 155                 TAIL <- "...\n\n"
 156             }
 157             cat("\nLabels:", paste(nms, collapse = " "), TAIL)
 158         } else {
 159             cat("1 Protein sequence in binary format stored in a vector.\n\n")
 160             cat("Sequence length:", nTot, "\n\n")
 161         }
 162     }
 163     if (nTot <= 1e6) {
 164         cat("Base composition:\n")
 165         print(round(base.freq(x), digits))
 166     } else cat("More than 1 million nucleotides: not printing base composition\n")
 167 }
 168
 169 as.proteinbin <- function(x, ...) UseMethod("as.proteinbin")
 170
 171 ._letters_ <-
 172     c("A","B","C","D","E","F","G","H","I","K","L","M","N","P",
 173       "Q","R","S","T","V","W","X","Y","Z","-")
 174 ._letters_alias_ <-
 175     list("?" = "X",
 176          "." = "-")
 177
 178 ### from http://en.wikipedia.org/wiki/Amino_acid and http://www.bioinformatics.org/sms2/iupac.html
 179 ._aa_3symbol_ <-
 180     list("A"="Ala",
 181          "B"="Asx",
 182          "R"="Arg",
 183          "N"="Asn",
 184          "D"="Asp",
 185          "C"="Cys",
 186          "E"="Glu",
 187          "Q"="Gln",
 188          "G"="Gly",
 189          "H"="His",
 190          "I"="Ile",
 191          "L"="Leu",
 192          "K"="Lys",
 193          "M"="Met",
 194          "F"="Phe",
 195          "P"="Pro",
 196          "S"="Ser",
 197          "T"="Thr",
 198          "W"="Trp",
 199          "Y"="Tyr",
 200          "V"="Val",
 201          "X"="Xaa",
 202          "Z"="Glx",
 203          "-"="Gap")
 204 ._aa_names_ <-
 205     list(
 206         "A"="Alanine",
 207         "B"="Apartic acid or Aspargine",
 208         "R"="Arginine",
 209         "N"="Asparagine",
 210         "D"="Aspartic acid",
 211         "C"="Cysteine",
 212         "E"="Glutamic acid",
 213         "Q"="Glutamine",
 214         "G"="Glycine",
 215         "H"="Histidine",
 216         "I"="Isoleucine",
 217         "L"="Leucine",
 218         "K"="Lysine",
 219         "M"="Methionine",
 220         "F"="Phenylalanine",
 221         "P"="Proline",
 222         "S"="Serine",
 223         "T"="Threonine",
 224         "W"="Tryptophan",
 225         "Y"="Tyrosine",
 226         "V"="Valine",
 227         "X"="Any amino acid",
 228         "Z"="Glutamine or Glutamic acid",
 229         "-"="Gap in alignment")
 230
 231 as.proteinbin.character <- function(x, ...)
 232 {
 233     ##
 234     if (any(nchar(x)>1)) {
 235         if (is.matrix(x)) {
 236             stop("Cannot convert a matrix of multiple characters into a sequence object")
 237         }
 238         x <- do.call(rbind,strsplit(x,""))
 239     }
 240     ans <- x
 241     for(alias in names(._letters_alias_)) {
 242         ans[ans==alias] <- ._letters_alias_[[alias]]
 243     }
 244     if (any(!(ans %in% ._letters_))) {
 245         stop("invalid characters for protein sequence")
 246     }
 247
 248     if (is.matrix(x)) {
 249         dim(ans) <- dim(x)
 250         dimnames(ans) <- dimnames(x)
 251     }
 252     class(ans) <- "proteinbin"
 253     ans
 254 }
 255
 256
 257 as.proteinbin.list <- function(x, ...)
 258 {
 259     obj <- lapply(x, as.proteinbin)
 260     class(obj) <- "proteinbin"
 261     obj
 262 }
 263
 264 ## they're already characters
 265 as.character.proteinbin <- function(x, ...)
 266 {
 267     f <- function(xx) {
 268         xx
 269     }
 270     if (is.list(x)) lapply(x, f) else f(x)
 271 }
 272
 273 base.freq <- function(x, freq = FALSE, all = FALSE)
 274 {
 275     f <- function(x) {
 276         ans <- table(x,deparse.level=0)
 277     }
 278
 279     if (is.list(x)) {
 280         BF <- rowSums(sapply(x, f))
 281         n <- sum(sapply(x, length))
 282     } else {
 283         n <- length(x)
 284         BF <- f(x)
 285     }
 286
 287     if (all) {
 288         if (!freq) BF <- BF / n
 289     } else {
 290         if (!freq) BF <- BF / sum(BF)
 291     }
 292     BF
 293 }
 294
 295
 296 ### where <- function(x, pattern)
 297 ### {
 298 ###     pat <- as.proteinbin(strsplit(pattern, NULL)[[1]])
 299 ###     p <- as.integer(length(pat))
 300 ###
 301 ###     foo <- function(x, pat, p) {
 302 ###         s <- as.integer(length(x))
 303 ###         if (s < p) stop("sequence shorter than the pattern")
 304 ###         ans <- .C("where", x, pat, s, p, integer(s), 0L,
 305 ###                   DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
 306 ###         n <- ans[[6]]
 307 ###         if (n) ans[[5]][seq_len(n)] - p + 2L else integer()
 308 ###     }
 309 ###
 310 ###     if (is.list(x)) return(lapply(x, foo, pat = pat, p = p))
 311 ###     if (is.matrix(x)) {
 312 ###         n <- nrow(x)
 313 ###         res <- vector("list", n)
 314 ###         for (i in seq_along(n))
 315 ###             res[[i]] <- foo(x[i, , drop = TRUE], pat, p)
 316 ###         names(res) <- rownames(x)
 317 ###         return(res)
 318 ###     }
 319 ###     foo(x, pat, p) # if x is a vector
 320 ### }
 321
 322