R/read.tree.R

   1 ## read.tree.R (2010-09-27)
   2
   3 ##   Read Tree Files in Parenthetic Format
   4
   5 ## Copyright 2002-2010 Emmanuel Paradis, Daniel Lawson and Klaus Schliep
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 tree.build <- function(tp)
  11 {
  12     add.internal <- function() {
  13         edge[j, 1] <<- current.node
  14         edge[j, 2] <<- current.node <<- node <<- node + 1L
  15         index[node] <<- j # set index
  16         j <<- j + 1L
  17     }
  18     add.terminal <- function() {
  19         edge[j, 1] <<- current.node
  20         edge[j, 2] <<- tip
  21         index[tip] <<- j # set index
  22         X <- unlist(strsplit(tpc[k], ":"))
  23         tip.label[tip] <<- X[1]
  24         edge.length[j] <<- as.numeric(X[2])
  25         k <<- k + 1L
  26         tip <<- tip + 1L
  27         j <<- j + 1L
  28     }
  29     go.down <- function() {
  30         l <- index[current.node]
  31         X <- unlist(strsplit(tpc[k], ":"))
  32         node.label[current.node - nb.tip] <<- X[1]
  33         edge.length[l] <<- as.numeric(X[2])
  34         k <<- k + 1L
  35         current.node <<- edge[l, 1]
  36     }
  37     if (!length(grep(",", tp))) {
  38         obj <- list(edge = matrix(c(2L, 1L), 1, 2))
  39         tp <- unlist(strsplit(tp, "[\\(\\):;]"))
  40         obj$edge.length <- as.numeric(tp[3])
  41         obj$Nnode <- 1L
  42         obj$tip.label <- tp[2]
  43         if (tp[4] != "") obj$node.label <- tp[4]
  44         class(obj) <- "phylo"
  45         return(obj)
  46     }
  47
  48     tpc <- unlist(strsplit(tp, "[\\(\\),;]"))
  49     tpc <- tpc[nzchar(tpc)]
  50     ## the following 2 lines are (slightly) faster than using gsub()
  51     tsp <- unlist(strsplit(tp, NULL))
  52     skeleton <- tsp[tsp %in% c("(", ")", ",", ";")]
  53     nsk <- length(skeleton)
  54     nb.node <- sum(skeleton == ")")
  55     nb.tip <- sum(skeleton == ",") + 1
  56     ## We will assume there is an edge at the root;
  57     ## if so, it will be removed and put into a vector
  58     nb.edge <- nb.node + nb.tip
  59     node.label <- character(nb.node)
  60     tip.label <- character(nb.tip)
  61
  62     edge.length <- numeric(nb.edge)
  63     edge <- matrix(0L, nb.edge, 2)
  64     current.node <- node <- as.integer(nb.tip + 1) # node number
  65     edge[nb.edge, 2] <- node
  66     index <- numeric(nb.edge + 1) # hash index to avoid which
  67     index[node] <- nb.edge
  68
  69     ## j: index of the line number of edge
  70     ## k: index of the line number of tpc
  71     ## tip: tip number
  72     j <- k <- tip <- 1L
  73
  74     for (i in 2:nsk) {
  75         if (skeleton[i] == "(") add.internal() # add an internal branch (on top)
  76         if (skeleton[i] == ",") {
  77             if (skeleton[i - 1] != ")") add.terminal() # add a terminal branch
  78         }
  79         if (skeleton[i] == ")") {
  80             if (skeleton[i - 1] == ",") { # add a terminal branch and go down one level
  81                 add.terminal()
  82                 go.down()
  83             }
  84             if (skeleton[i - 1] == ")") go.down() # go down one level
  85         }
  86     }
  87
  88     edge <- edge[-nb.edge, ]
  89     obj <- list(edge = edge, Nnode = nb.node, tip.label = tip.label)
  90     root.edge <- edge.length[nb.edge]
  91     edge.length <- edge.length[-nb.edge]
  92     if (!all(is.na(edge.length))) # added 2005-08-18
  93         obj$edge.length <- edge.length
  94     if (is.na(node.label[1])) node.label[1] <- ""
  95     if (any(nzchar(node.label))) obj$node.label <- node.label
  96     if (!is.na(root.edge)) obj$root.edge <- root.edge
  97     class(obj) <- "phylo"
  98     obj
  99 }
 100
 101 read.tree <- function(file = "", text = NULL, tree.names = NULL, skip = 0,
 102     comment.char = "#", keep.multi = FALSE, ...)
 103 {
 104     unname <- function(treetext) {
 105         nc <- nchar(treetext)
 106         tstart <- 1
 107         while (substr(treetext, tstart, tstart) != "(" && tstart <= nc)
 108             tstart <- tstart + 1
 109         if (tstart > 1)
 110             return(c(substr(treetext, 1, tstart - 1),
 111                      substr(treetext, tstart, nc)))
 112         return(c("", treetext))
 113     }
 114     if (!is.null(text)) {
 115         if (!is.character(text))
 116           stop("argument `text' must be of mode character")
 117         tree <- text
 118     } else {
 119         tree <- scan(file = file, what = "", sep = "\n", quiet = TRUE,
 120                      skip = skip, comment.char = comment.char, ...)
 121     }
 122     ## Suggestion from Eric Durand and Nicolas Bortolussi (added 2005-08-17):
 123     if (identical(tree, character(0))) {
 124         warning("empty character string.")
 125         return(NULL)
 126     }
 127     tree <- gsub("[ \t]", "", tree)
 128     tree <- unlist(strsplit(tree, NULL))
 129     y <- which(tree == ";")
 130     Ntree <- length(y)
 131     x <- c(1, y[-Ntree] + 1)
 132     ## Suggestion from Olivier François (added 2006-07-15):
 133     if (is.na(y[1])) return(NULL)
 134     STRING <- character(Ntree)
 135     for (i in 1:Ntree)
 136         STRING[i] <- paste(tree[x[i]:y[i]], sep = "", collapse = "")
 137
 138     tmp <- unlist(lapply(STRING, unname))
 139     tmpnames <- tmp[c(TRUE, FALSE)]
 140     STRING <- tmp[c(FALSE, TRUE)]
 141     if (is.null(tree.names) && any(nzchar(tmpnames)))
 142         tree.names <- tmpnames
 143
 144     colon <- grep(":", STRING)
 145     if (!length(colon)) {
 146         obj <- lapply(STRING, clado.build)
 147     } else if (length(colon) == Ntree) {
 148         obj <- lapply(STRING, tree.build)
 149     } else {
 150         obj <- vector("list", Ntree)
 151         obj[colon] <- lapply(STRING[colon], tree.build)
 152         nocolon <- (1:Ntree)[!1:Ntree %in% colon]
 153         obj[nocolon] <- lapply(STRING[nocolon], clado.build)
 154     }
 155     for (i in 1:Ntree) {
 156         ## Check here that the root edge is not incorrectly represented
 157         ## in the object of class "phylo" by simply checking that there
 158         ## is a bifurcation at the root
 159         ROOT <- length(obj[[i]]$tip.label) + 1
 160         if(sum(obj[[i]]$edge[, 1] == ROOT) == 1 && dim(obj[[i]]$edge)[1] > 1)
 161             stop(paste("There is apparently two root edges in your file: cannot read tree file.\n  Reading Newick file aborted at tree no.", i, sep = ""))
 162     }
 163     if (Ntree == 1 && !keep.multi) obj <- obj[[1]] else {
 164         if (!is.null(tree.names)) names(obj) <- tree.names
 165         class(obj) <- "multiPhylo"
 166     }
 167     obj
 168 }