R/read.dna.R

   1 ## read.dna.R (2012-05-03)
   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.dna <- function(file, format = "interleaved", skip = 0,
  11                      nlines = 0, comment.char = "#",
  12                      as.character = FALSE, as.matrix = NULL)
  13 {
  14     findFirstNucleotide <- function(x) {
  15         ## actually find the 1st non-blank character
  16         ## just in case: pat.base <- "[-AaCcGgTtUuMmRrWwSsYyKkVvHhDdBbNn?]{10}"
  17         tmp <- regexpr("[[:blank:]]+", x[1]) # consider only a single string
  18         tmp[1] + attr(tmp, "match.length")
  19     }
  20     getTaxaNames <- function(x) {
  21         x <- sub("^['\" ]+", "", x) # remove the leading quotes and spaces
  22         x <- sub("['\" ]+$", "", x) #   "     "  trailing  "     "    "
  23         x
  24     }
  25     getNucleotide <- function(x) {
  26         x <- gsub(" ", "", x)
  27         x <- strsplit(x, NULL)
  28         tolower(unlist(x))
  29     }
  30     formats <- c("interleaved", "sequential", "fasta", "clustal")
  31     format <- match.arg(format, formats)
  32     phylip <- if (format %in% formats[1:2]) TRUE else FALSE
  33     X <- scan(file = file, what = "", sep = "\n", quiet = TRUE,
  34               skip = skip, nlines = nlines, comment.char = comment.char)
  35
  36     if (phylip) {
  37         ## need to remove the possible leading spaces and/or tabs in the first line
  38         fl <- gsub("^[[:blank:]]+", "", X[1])
  39         fl <- as.numeric(unlist(strsplit(fl, "[[:blank:]]+")))
  40         if (length(fl) != 2 || any(is.na(fl)))
  41             stop("the first line of the file must contain the dimensions of the data")
  42         n <- fl[1]
  43         s <- fl[2]
  44         obj <- matrix("", n, s)
  45         X <- X[-1]
  46     }
  47     switch(format,
  48     "interleaved" = {
  49         start.seq <- findFirstNucleotide(X[1])
  50         one2n <- 1:n
  51         taxa <- getTaxaNames(substr(X[one2n], 1, start.seq - 1))
  52         X[one2n] <- substr(X[one2n], start.seq, nchar(X[one2n]))
  53         nl <- length(X)
  54         for (i in one2n)
  55             obj[i, ] <- getNucleotide(X[seq(i, nl, n)])
  56     },
  57     "sequential" = {
  58         taxa <- character(n)
  59         j <- 1L # line number
  60         for (i in 1:n) {
  61             start.seq <- findFirstNucleotide(X[j])
  62             taxa[i] <- getTaxaNames(substr(X[j], 1, start.seq - 1))
  63             sequ <- getNucleotide(substr(X[j], start.seq, nchar(X[j])))
  64             j <- j + 1L
  65             while (length(sequ) < s) {
  66                 sequ <- c(sequ, getNucleotide(X[j]))
  67                 j <- j + 1L
  68             }
  69             obj[i, ] <- sequ
  70         }
  71         taxa <- getTaxaNames(taxa)
  72     },
  73     "fasta" = {
  74         start <- grep("^ {0,}>", X)
  75         taxa <- X[start]
  76         taxa <- sub("^ {0,}>", "", taxa) # remove the hook and the spaces before
  77         taxa <- getTaxaNames(taxa)
  78         n <- length(taxa)
  79         obj <- vector("list", n)
  80         start <- c(start, length(X) + 1) # this avoids the following to crash when `i = n'
  81         for (i in 1:n)
  82             obj[[i]] <- getNucleotide(X[(start[i] + 1):(start[i + 1] - 1)])
  83     },
  84     "clustal" = {
  85         X <- X[-1] # drop the line with "Clustal bla bla..."
  86         ## find where the 1st sequence starts
  87         start.seq <- findFirstNucleotide(X[1])
  88         ## find the lines with *********....
  89         nspaces <- paste("^ {", start.seq - 1, "}", sep = "", collapse = "")
  90         stars <- grep(nspaces, X)
  91         ## we now know how many sequences in the file:
  92         n <- stars[1] - 1
  93         taxa <- getTaxaNames(substr(X[1:n], 1, start.seq - 1))
  94         ## need to remove the sequence names before getting the sequences:
  95         X <- substr(X, start.seq, nchar(X))
  96         nl <- length(X)
  97         ## find the length of the 1st sequence:
  98         tmp <- getNucleotide(X[seq(1, nl, n + 1)])
  99         s <- length(tmp)
 100         obj <- matrix("", n, s)
 101         obj[1, ] <- tmp
 102         for (i in 2:n)
 103             obj[i, ] <- getNucleotide(X[seq(i, nl, n + 1)])
 104     })
 105
 106     if (format != "fasta") {
 107         rownames(obj) <- taxa
 108     } else {
 109         names(obj) <- taxa
 110         LENGTHS <- unique(unlist(lapply(obj, length)))
 111         allSameLength <- length(LENGTHS) == 1
 112         if (is.logical(as.matrix)) {
 113             if (as.matrix && !allSameLength)
 114                 stop("sequences in FASTA file not of the same length")
 115         } else {
 116             as.matrix <- allSameLength
 117         }
 118         if (as.matrix) {
 119             obj <- matrix(unlist(obj), ncol = LENGTHS, byrow = TRUE)
 120             rownames(obj) <- taxa
 121         }
 122     }
 123     if (!as.character) obj <- as.DNAbin(obj)
 124     obj
 125 }