fix in read.dna and change in write.dna

[ape.git] / R / read.dna.R

diff --git a/R/read.dna.R b/R/read.dna.R
index 3e2b980670d872529f78eae740215d5d1582c279..ca4da7a4daac444adbfab40c413d6ff08fa5f72d 100644 (file)
--- a/R/read.dna.R
+++ b/R/read.dna.R
@@ -1,4 +1,4 @@
-## read.dna.R (2012-04-26)
+## read.dna.R (2012-05-03)
 
 ##   Read DNA Sequences in a File
 
@@ -8,9 +8,15 @@
 ## See the file ../COPYING for licensing issues.
 
 read.dna <- function(file, format = "interleaved", skip = 0,
-                     nlines = 0, comment.char = "#", seq.names = NULL,
+                     nlines = 0, comment.char = "#",
                      as.character = FALSE, as.matrix = NULL)
 {
+    findFirstNucleotide <- function(x) {
+        ## actually find the 1st non-blank character
+        ## just in case: pat.base <- "[-AaCcGgTtUuMmRrWwSsYyKkVvHhDdBbNn?]{10}"
+        tmp <- regexpr("[[:blank:]]+", x[1]) # consider only a single string
+        tmp[1] + attr(tmp, "match.length")
+    }
     getTaxaNames <- function(x) {
         x <- sub("^['\" ]+", "", x) # remove the leading quotes and spaces
         x <- sub("['\" ]+$", "", x) #   "     "  trailing  "     "    "
@@ -19,14 +25,14 @@ read.dna <- function(file, format = "interleaved", skip = 0,
     getNucleotide <- function(x) {
         x <- gsub(" ", "", x)
         x <- strsplit(x, NULL)
-        unlist(x)
+        tolower(unlist(x))
     }
     formats <- c("interleaved", "sequential", "fasta", "clustal")
     format <- match.arg(format, formats)
     phylip <- if (format %in% formats[1:2]) TRUE else FALSE
     X <- scan(file = file, what = "", sep = "\n", quiet = TRUE,
               skip = skip, nlines = nlines, comment.char = comment.char)
-    pat.base <- "[-AaCcGgTtUuMmRrWwSsYyKkVvHhDdBbNn?]{10}"
+
     if (phylip) {
         ## need to remove the possible leading spaces and/or tabs in the first line
         fl <- gsub("^[[:blank:]]+", "", X[1])
@@ -38,56 +44,53 @@ read.dna <- function(file, format = "interleaved", skip = 0,
         obj <- matrix("", n, s)
         X <- X[-1]
     }
-    if (format == "interleaved") {
-        start.seq <- regexpr(pat.base, X[1])[1]
-        if (is.null(seq.names))
-            seq.names <- getTaxaNames(substr(X[1:n], 1, start.seq - 1))
-        X[1:n] <- substr(X[1:n], start.seq, nchar(X[1:n]))
+    switch(format,
+    "interleaved" = {
+        start.seq <- findFirstNucleotide(X[1])
+        one2n <- 1:n
+        taxa <- getTaxaNames(substr(X[one2n], 1, start.seq - 1))
+        X[one2n] <- substr(X[one2n], start.seq, nchar(X[one2n]))
         nl <- length(X)
-        for (i in 1:n)
+        for (i in one2n)
             obj[i, ] <- getNucleotide(X[seq(i, nl, n)])
-    }
-    if (format == "sequential") {
+    },
+    "sequential" = {
         taxa <- character(n)
-        j <- 1
+        j <- 1L # line number
         for (i in 1:n) {
-            start.seq <- regexpr(pat.base, X[j])[1]
-            taxa[i] <- substr(X[j], 1, start.seq - 1)
+            start.seq <- findFirstNucleotide(X[j])
+            taxa[i] <- getTaxaNames(substr(X[j], 1, start.seq - 1))
             sequ <- getNucleotide(substr(X[j], start.seq, nchar(X[j])))
-            j <- j + 1
+            j <- j + 1L
             while (length(sequ) < s) {
                 sequ <- c(sequ, getNucleotide(X[j]))
-                j <- j + 1
+                j <- j + 1L
             }
             obj[i, ] <- sequ
         }
-        if (is.null(seq.names)) seq.names <- getTaxaNames(taxa)
-    }
-    if (format == "fasta") {
+        taxa <- getTaxaNames(taxa)
+    },
+    "fasta" = {
         start <- grep("^ {0,}>", X)
         taxa <- X[start]
+        taxa <- sub("^ {0,}>", "", taxa) # remove the hook and the spaces before
+        taxa <- getTaxaNames(taxa)
         n <- length(taxa)
         obj <- vector("list", n)
-        if (is.null(seq.names)) {
-            taxa <- sub("^ {0,}>", "", taxa) # remove the hook and the spaces before
-            seq.names <- getTaxaNames(taxa)
-        }
         start <- c(start, length(X) + 1) # this avoids the following to crash when `i = n'
         for (i in 1:n)
             obj[[i]] <- getNucleotide(X[(start[i] + 1):(start[i + 1] - 1)])
-    }
-    if (format == "clustal") {
-        X <- X[-1]
+    },
+    "clustal" = {
+        X <- X[-1] # drop the line with "Clustal bla bla..."
         ## find where the 1st sequence starts
-        start.seq <- regexpr(pat.base, X[1])[1]
+        start.seq <- findFirstNucleotide(X[1])
         ## find the lines with *********....
         nspaces <- paste("^ {", start.seq - 1, "}", sep = "", collapse = "")
         stars <- grep(nspaces, X)
         ## we now know how many sequences in the file:
         n <- stars[1] - 1
-        ## get the sequence names in the same way than "interleaved":
-        if (is.null(seq.names))
-            seq.names <- getTaxaNames(substr(X[1:n], 1, start.seq - 1))
+        taxa <- getTaxaNames(substr(X[1:n], 1, start.seq - 1))
         ## need to remove the sequence names before getting the sequences:
         X <- substr(X, start.seq, nchar(X))
         nl <- length(X)
@@ -98,13 +101,12 @@ read.dna <- function(file, format = "interleaved", skip = 0,
         obj[1, ] <- tmp
         for (i in 2:n)
             obj[i, ] <- getNucleotide(X[seq(i, nl, n + 1)])
-    }
+    })
+
     if (format != "fasta") {
-        rownames(obj) <- seq.names
-        obj <- tolower(obj)
+        rownames(obj) <- taxa
     } else {
-        names(obj) <- seq.names
-        obj <- lapply(obj, tolower)
+        names(obj) <- taxa
         LENGTHS <- unique(unlist(lapply(obj, length)))
         allSameLength <- length(LENGTHS) == 1
         if (is.logical(as.matrix)) {
@@ -115,7 +117,7 @@ read.dna <- function(file, format = "interleaved", skip = 0,
         }
         if (as.matrix) {
             obj <- matrix(unlist(obj), ncol = LENGTHS, byrow = TRUE)
-            rownames(obj) <- seq.names
+            rownames(obj) <- taxa
         }
     }
     if (!as.character) obj <- as.DNAbin(obj)