-## read.dna.R (2010-05-17)
+## read.dna.R (2013-04-02)
## Read DNA Sequences in a File
-## Copyright 2003-2010 Emmanuel Paradis
+## Copyright 2003-2013 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
+read.FASTA <- function(file)
+{
+ sz <- file.info(file)$size
+ x <- readBin(file, "raw", sz)
+ icr <- which(x == as.raw(0x0d)) # CR
+ if (length(icr)) x <- x[-icr]
+ res <- .Call("rawStreamToDNAbin", x, PACKAGE = "ape")
+ names(res) <- sub("^ +", "", names(res)) # to permit phylosim
+ class(res) <- "DNAbin"
+ res
+}
+
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 " " "
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 in the first line
- fl <- gsub("^ +", "", X[1])
- fl <- as.numeric(unlist(strsplit(fl, " +")))
- if (length(fl) != 2 || any(is.na(fl)))
- stop("the first line of the file must contain the dimensions of the data")
- n <- fl[1]
- s <- fl[2]
- 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]))
- nl <- length(X)
- for (i in 1:n)
- obj[i, ] <- getNucleotide(X[seq(i, nl, n)])
- }
- if (format == "sequential") {
- taxa <- character(n)
- j <- 1
- for (i in 1:n) {
- start.seq <- regexpr(pat.base, X[j])[1]
- taxa[i] <- substr(X[j], 1, start.seq - 1)
- sequ <- getNucleotide(substr(X[j], start.seq, nchar(X[j])))
- j <- j + 1
- while (length(sequ) < s) {
- sequ <- c(sequ, getNucleotide(X[j]))
- j <- j + 1
- }
- obj[i, ] <- sequ
- }
- if (is.null(seq.names)) seq.names <- getTaxaNames(taxa)
- }
if (format == "fasta") {
- start <- grep("^ {0,}>", X)
- taxa <- X[start]
- 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)
+ obj <- read.FASTA(file)
+ } else {
+ X <- scan(file = file, what = "", sep = "\n", quiet = TRUE,
+ skip = skip, nlines = nlines, comment.char = comment.char)
+ if (format %in% formats[1:2]) {
+ ## need to remove the possible leading spaces and/or tabs in the first line
+ fl <- gsub("^[[:blank:]]+", "", X[1])
+ fl <- as.numeric(unlist(strsplit(fl, "[[:blank:]]+")))
+ if (length(fl) != 2 || any(is.na(fl)))
+ stop("the first line of the file must contain the dimensions of the data")
+ n <- fl[1]
+ s <- fl[2]
+ obj <- matrix("", n, s)
+ X <- X[-1]
}
- 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]
- ## find where the 1st sequence starts
- start.seq <- regexpr(pat.base, X[1])[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))
- ## need to remove the sequence names before getting the sequences:
- X <- substr(X, start.seq, nchar(X))
- nl <- length(X)
- ## find the length of the 1st sequence:
- tmp <- getNucleotide(X[seq(1, nl, n + 1)])
- s <- length(tmp)
- obj <- matrix("", n, s)
- obj[1, ] <- tmp
- for (i in 2:n)
- obj[i, ] <- getNucleotide(X[seq(i, nl, n + 1)])
+ 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 one2n)
+ obj[i, ] <- getNucleotide(X[seq(i, nl, n)])
+ },
+ "sequential" = {
+ taxa <- character(n)
+ j <- 1L # line number
+ for (i in 1:n) {
+ 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 + 1L
+ while (length(sequ) < s) {
+ sequ <- c(sequ, getNucleotide(X[j]))
+ j <- j + 1L
+ }
+ obj[i, ] <- sequ
+ }
+ taxa <- getTaxaNames(taxa)
+ },
+ "clustal" = {
+ X <- X[-1] # drop the line with "Clustal bla bla..."
+ ## find where the 1st sequence starts
+ 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
+ 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)
+ ## find the length of the 1st sequence:
+ tmp <- getNucleotide(X[seq(1, nl, n + 1)])
+ s <- length(tmp)
+ obj <- matrix("", n, s)
+ 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
+ if (!as.character) obj <- as.DNAbin(obj)
} else {
- names(obj) <- seq.names
- obj <- lapply(obj, tolower)
LENGTHS <- unique(unlist(lapply(obj, length)))
allSameLength <- length(LENGTHS) == 1
- if (is.logical(as.matrix) && as.matrix && !allSameLength)
- stop("sequences in FASTA file not of the same length")
- if (is.null(as.matrix) && allSameLength)
- as.matrix <- TRUE
+ if (is.logical(as.matrix)) {
+ if (as.matrix && !allSameLength)
+ stop("sequences in FASTA file not of the same length")
+ } else {
+ as.matrix <- allSameLength
+ }
if (as.matrix) {
- obj <- matrix(unlist(obj), ncol = LENGTHS, byrow = TRUE)
- rownames(obj) <- seq.names
+ taxa <- names(obj)
+ n <- length(obj)
+ y <- matrix(as.raw(0), n, LENGTHS)
+ for (i in seq_len(n)) y[i, ] <- obj[[i]]
+ obj <- y
+ rownames(obj) <- taxa
+ class(obj) <- "DNAbin"
}
+ if (as.character) obj <- as.character(obj)
}
- if (!as.character) obj <- as.DNAbin(obj)
obj
}
projects / ape.git / blobdiff