R/read.dna.R

   1 ## read.dna.R (2013-01-04)
   2
   3 ##   Read DNA Sequences in a File
   4
   5 ## Copyright 2003-2012 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 read.FASTA <- function(file)
  11 {
  12     sz <- file.info(file)$size
  13     x <- readBin(file, "raw", sz)
  14     if (Sys.info()[1] == "Windows") {
  15         icr <- which(x == as.raw(0x0d)) # CR
  16         x <- x[-icr]
  17     }
  18     res <- .Call("rawStreamToDNAbin", x, PACKAGE = "ape")
  19     names(res) <- sub("^ +", "", names(res)) # to permit phylosim
  20     class(res) <- "DNAbin"
  21     res
  22 }
  23
  24 read.dna <- function(file, format = "interleaved", skip = 0,
  25                      nlines = 0, comment.char = "#",
  26                      as.character = FALSE, as.matrix = NULL)
  27 {
  28     findFirstNucleotide <- function(x) {
  29         ## actually find the 1st non-blank character
  30         ## just in case: pat.base <- "[-AaCcGgTtUuMmRrWwSsYyKkVvHhDdBbNn?]{10}"
  31         tmp <- regexpr("[[:blank:]]+", x[1]) # consider only a single string
  32         tmp[1] + attr(tmp, "match.length")
  33     }
  34     getTaxaNames <- function(x) {
  35         x <- sub("^['\" ]+", "", x) # remove the leading quotes and spaces
  36         x <- sub("['\" ]+$", "", x) #   "     "  trailing  "     "    "
  37         x
  38     }
  39     getNucleotide <- function(x) {
  40         x <- gsub(" ", "", x)
  41         x <- strsplit(x, NULL)
  42         tolower(unlist(x))
  43     }
  44     formats <- c("interleaved", "sequential", "fasta", "clustal")
  45     format <- match.arg(format, formats)
  46     if (format == "fasta") {
  47         obj <- read.FASTA(file)
  48     } else {
  49         X <- scan(file = file, what = "", sep = "\n", quiet = TRUE,
  50                   skip = skip, nlines = nlines, comment.char = comment.char)
  51         if (format %in% formats[1:2]) {
  52             ## need to remove the possible leading spaces and/or tabs in the first line
  53             fl <- gsub("^[[:blank:]]+", "", X[1])
  54             fl <- as.numeric(unlist(strsplit(fl, "[[:blank:]]+")))
  55             if (length(fl) != 2 || any(is.na(fl)))
  56                 stop("the first line of the file must contain the dimensions of the data")
  57             n <- fl[1]
  58             s <- fl[2]
  59             obj <- matrix("", n, s)
  60             X <- X[-1]
  61         }
  62         switch(format,
  63                "interleaved" = {
  64                    start.seq <- findFirstNucleotide(X[1])
  65                    one2n <- 1:n
  66                    taxa <- getTaxaNames(substr(X[one2n], 1, start.seq - 1))
  67                    X[one2n] <- substr(X[one2n], start.seq, nchar(X[one2n]))
  68                    nl <- length(X)
  69                    for (i in one2n)
  70                        obj[i, ] <- getNucleotide(X[seq(i, nl, n)])
  71                },
  72                "sequential" = {
  73                    taxa <- character(n)
  74                    j <- 1L # line number
  75                    for (i in 1:n) {
  76                        start.seq <- findFirstNucleotide(X[j])
  77                        taxa[i] <- getTaxaNames(substr(X[j], 1, start.seq - 1))
  78                        sequ <- getNucleotide(substr(X[j], start.seq, nchar(X[j])))
  79                        j <- j + 1L
  80                        while (length(sequ) < s) {
  81                            sequ <- c(sequ, getNucleotide(X[j]))
  82                            j <- j + 1L
  83                        }
  84                        obj[i, ] <- sequ
  85                    }
  86                    taxa <- getTaxaNames(taxa)
  87                },
  88                "clustal" = {
  89                    X <- X[-1] # drop the line with "Clustal bla bla..."
  90                    ## find where the 1st sequence starts
  91                    start.seq <- findFirstNucleotide(X[1])
  92                    ## find the lines with *********....
  93                    nspaces <- paste("^ {", start.seq - 1, "}", sep = "", collapse = "")
  94                    stars <- grep(nspaces, X)
  95                    ## we now know how many sequences in the file:
  96                    n <- stars[1] - 1
  97                    taxa <- getTaxaNames(substr(X[1:n], 1, start.seq - 1))
  98                    ## need to remove the sequence names before getting the sequences:
  99                    X <- substr(X, start.seq, nchar(X))
 100                    nl <- length(X)
 101                    ## find the length of the 1st sequence:
 102                    tmp <- getNucleotide(X[seq(1, nl, n + 1)])
 103                    s <- length(tmp)
 104                    obj <- matrix("", n, s)
 105                    obj[1, ] <- tmp
 106                    for (i in 2:n)
 107                        obj[i, ] <- getNucleotide(X[seq(i, nl, n + 1)])
 108                })
 109     }
 110     if (format != "fasta") {
 111         rownames(obj) <- taxa
 112         if (!as.character) obj <- as.DNAbin(obj)
 113     } else {
 114         LENGTHS <- unique(unlist(lapply(obj, length)))
 115         allSameLength <- length(LENGTHS) == 1
 116         if (is.logical(as.matrix)) {
 117             if (as.matrix && !allSameLength)
 118                 stop("sequences in FASTA file not of the same length")
 119         } else {
 120             as.matrix <- allSameLength
 121         }
 122         if (as.matrix) {
 123             taxa <- names(obj)
 124             n <- length(obj)
 125             y <- matrix(as.raw(0), n, LENGTHS)
 126             for (i in seq_len(n)) y[i, ] <- obj[[i]]
 127             obj <- y
 128             rownames(obj) <- taxa
 129             class(obj) <- "DNAbin"
 130         }
 131         if (as.character) obj <- as.character(obj)
 132     }
 133     obj
 134 }