X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Fread.dna.R;h=ab8051e023afeb2c48c4613348fe360c4dce224c;hb=b0d1251527d8dd48ca1703b1cfaf217f413eda0e;hp=98ebd65ca71c67e2812a50dec6bf118e3b023f5e;hpb=dfabbb4721bc499ee9e5ee9aadcfd59ff3d4f223;p=ape.git diff --git a/R/read.dna.R b/R/read.dna.R index 98ebd65..ab8051e 100644 --- a/R/read.dna.R +++ b/R/read.dna.R @@ -1,16 +1,34 @@ -## 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 " " " @@ -19,103 +37,96 @@ 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 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 }