R/read.dna.R

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