R/read.nexus.R

   1 ## read.nexus.R (2012-09-28)
   2
   3 ##   Read Tree File in Nexus Format
   4
   5 ## Copyright 2003-2012 Emmanuel Paradis and 2010 Klaus Schliep
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 .treeBuildWithTokens <- function(x)
  11 {
  12     phy <- .Call("treeBuildWithTokens", x, PACKAGE = "ape")
  13     dim(phy[[1]]) <- c(length(phy[[1]])/2, 2)
  14     nms <- c("edge", "edge.length", "Nnode", "node.label", "root.edge")
  15     if (length(phy) == 4) nms <- nms[-5]
  16     names(phy) <- nms
  17     if (all(phy$node.label == "")) phy$node.label <- NULL
  18     class(phy) <- "phylo"
  19     attr(phy, "order") <- "cladewise"
  20     phy
  21 }
  22
  23 clado.build <- function(tp)
  24 {
  25     add.internal <- function() {
  26         edge[j, 1] <<- current.node
  27         node <<- node + 1
  28         edge[j, 2] <<- current.node <<- node
  29         index[node] <<- j # set index
  30         j <<- j + 1
  31     }
  32     add.terminal <- function() {
  33         edge[j, 1] <<- current.node
  34         edge[j, 2] <<- tip
  35         index[tip] <<- j # set index
  36         tip.label[tip] <<- tpc[k]
  37         k <<- k + 1
  38         tip <<- tip + 1
  39         j <<- j + 1
  40     }
  41     go.down <- function() {
  42         l <- index[current.node]
  43         node.label[current.node - nb.tip] <<- tpc[k]
  44         k <<- k + 1
  45         current.node <<- edge[l, 1]
  46     }
  47     if (!length(grep(",", tp))) {
  48         obj <- list(edge = matrix(c(2, 1), 1, 2), Nnode = 1)
  49         tp <- unlist(strsplit(tp, "[\\(\\);]"))
  50         obj$tip.label <- tp[2]
  51         if (tp[3] != "") obj$node.label <- tp[3]
  52         class(obj) <- "phylo"
  53         return(obj)
  54     }
  55     tsp <- unlist(strsplit(tp, NULL))
  56     tp <- gsub(")", ")NA", tp)
  57     tp <- gsub(" ", "", tp)
  58     tpc <- unlist(strsplit(tp, "[\\(\\),;]"))
  59     tpc <- tpc[tpc != ""]
  60     skeleton <- tsp[tsp == "(" | tsp == ")" | tsp == "," | tsp == ";"]
  61     nsk <- length(skeleton)
  62     nb.node <- length(skeleton[skeleton == ")"])
  63     nb.tip <- length(skeleton[skeleton == ","]) + 1
  64     ## We will assume there is an edge at the root;
  65     ## if so, it will be removed and put in a vector
  66     nb.edge <- nb.node + nb.tip
  67     node.label <- character(nb.node)
  68     tip.label <- character(nb.tip)
  69
  70     edge <- matrix(NA, nb.edge, 2)
  71     current.node <- node <- nb.tip + 1 # node number
  72     edge[nb.edge, 1] <- 0    # see comment above
  73     edge[nb.edge, 2] <- node #
  74
  75     index <- numeric(nb.edge + 1)
  76     index[node] <- nb.edge
  77     ## j: index of the line number of edge
  78     ## k: index of the line number of tpc
  79     ## tip: tip number
  80     j <- k <- tip <- 1
  81     for (i in 2:nsk) {
  82         if (skeleton[i] == "(") add.internal()      # add an internal branch (on top)
  83         if (skeleton[i] == ",") {
  84             if (skeleton[i - 1] != ")") add.terminal()   # add a terminal branch
  85         }
  86         if (skeleton[i] == ")") {
  87             if (skeleton[i - 1] == ",") {   # add a terminal branch and go down one level
  88                 add.terminal()
  89                 go.down()
  90             }
  91             if (skeleton[i - 1] == ")") go.down()   # go down one level
  92         }
  93     }
  94     edge <- edge[-nb.edge, ]
  95     obj <- list(edge = edge, tip.label = tip.label,
  96                 Nnode = nb.node, node.label = node.label)
  97     obj$node.label <-
  98         if (all(obj$node.label == "NA")) NULL
  99         else gsub("^NA", "", obj$node.label)
 100     class(obj) <- "phylo"
 101     attr(obj, "order") <- "cladewise"
 102     obj
 103 }
 104
 105 read.nexus <- function(file, tree.names = NULL)
 106 {
 107     X <- scan(file = file, what = "", sep = "\n", quiet = TRUE)
 108     ## remove all comments
 109     ## (this might not work if there are square brackets within the comments)
 110     LEFT <- grep("\\[", X)
 111     RIGHT <- grep("\\]", X)
 112     if (length(LEFT)) { # in case there are no comments at all
 113         w <- LEFT == RIGHT
 114         if (any(w)) { # in case all comments use at least 2 lines
 115             s <- LEFT[w]
 116             X[s] <- gsub("\\[[^]]*\\]", "", X[s])
 117             ## The above regexp was quite tough to find: it makes
 118             ## possible to delete series of comments on the same line:
 119             ##       ...[...]xxx[...]...
 120             ## without deleting the "xxx". This regexp is in three parts:
 121             ##       \\[      [^]]*       \\]
 122             ## where [^]]* means "any character, except "]", repeated zero
 123             ## or more times" (note that the ']' is not escaped here).
 124             ## The previous version was:
 125             ##       X[s] <- gsub("\\[.*\\]", "", X[s])
 126             ## which deleted the "xxx". (EP  2008-06-24)
 127         }
 128         w <- !w
 129         if (any(w)) {
 130             s <- LEFT[w]
 131             X[s] <- gsub("\\[.*", "", X[s])
 132             sb <- RIGHT[w]
 133             X[sb] <- gsub(".*\\]", "", X[sb])
 134             if (any(s < sb - 1))
 135                 X <- X[-unlist(mapply(":", (s + 1), (sb - 1)))]
 136         }
 137     }
 138     endblock <- grep("END;|ENDBLOCK;", X, ignore.case = TRUE)
 139     semico <- grep(";", X)
 140     i1 <- grep("BEGIN TREES;", X, ignore.case = TRUE)
 141     i2 <- grep("TRANSLATE", X, ignore.case = TRUE)
 142     translation <- if (length(i2) == 1 && i2 > i1) TRUE else FALSE
 143     if (translation) {
 144         end <- semico[semico > i2][1]
 145         x <- X[(i2 + 1):end] # assumes there's a 'new line' after "TRANSLATE"
 146         ## x <- gsub("TRANSLATE", "", x, ignore.case = TRUE)
 147         x <- unlist(strsplit(x, "[,; \t]"))
 148         x <- x[nzchar(x)]
 149         TRANS <- matrix(x, ncol = 2, byrow = TRUE)
 150         TRANS[, 2] <- gsub("['\"]", "", TRANS[, 2])
 151         n <- dim(TRANS)[1]
 152     }
 153     start <-
 154         if (translation) semico[semico > i2][1] + 1
 155         else semico[semico > i1][1]
 156     end <- endblock[endblock > i1][1] - 1
 157     tree <- X[start:end]
 158     rm(X)
 159 ###    tree <- gsub("^.*= *", "", tree)
 160     ## check whether there are empty lines from the above manips:
 161     tree <- tree[tree != ""]
 162     semico <- grep(";", tree)
 163     Ntree <- length(semico) # provisional -- some ";" may actually mark end of commands
 164     ## are some trees on several lines?
 165     ## -- this actually 'packs' all characters ending with a ";" in a single string --
 166     if (Ntree == 1 && length(tree) > 1) STRING <- paste(tree, collapse = "") else {
 167         if (any(diff(semico) != 1)) {
 168             STRING <- character(Ntree)
 169             s <- c(1, semico[-Ntree] + 1)
 170             j <- mapply(":", s, semico)
 171             if (is.list(j)) {
 172                 for (i in 1:Ntree)
 173                     STRING[i] <- paste(tree[j[[i]]], collapse = "")
 174             } else {
 175                 for (i in 1:Ntree)
 176                     STRING[i] <- paste(tree[j[, i]], collapse = "")
 177             }
 178         } else STRING <- tree
 179     }
 180     rm(tree)
 181     ## exclude the possible command lines ending with ";":
 182     STRING <- STRING[grep("^[[:blank:]]*tree.*= *", STRING, ignore.case = TRUE)]
 183     Ntree <- length(STRING) # update Ntree
 184     ## get the tree names:
 185     nms.trees <- sub(" *= *.*", "", STRING) # only the first occurence of "="
 186     nms.trees <- sub("^ *tree *", "", nms.trees, ignore.case = TRUE)
 187     STRING <- sub("^.*= *", "", STRING) # delete title and 'TREE' command with 'sub'
 188     STRING <- gsub(" ", "", STRING) # delete all white spaces
 189     colon <- grep(":", STRING)
 190     if (!length(colon)) {
 191         trees <- lapply(STRING, clado.build)
 192     } else if (length(colon) == Ntree) {
 193         trees <-
 194             if (translation) lapply(STRING, .treeBuildWithTokens)
 195             else lapply(STRING, tree.build)
 196     } else {
 197         trees <- vector("list", Ntree)
 198         trees[colon] <- lapply(STRING[colon], tree.build)
 199         nocolon <- (1:Ntree)[!1:Ntree %in% colon]
 200         trees[nocolon] <- lapply(STRING[nocolon], clado.build)
 201         if (translation) {
 202             for (i in 1:Ntree) {
 203                 tr <- trees[[i]]
 204                 for (j in 1:n) {
 205                     ind <- which(tr$tip.label[j] == TRANS[, 1])
 206                     tr$tip.label[j] <- TRANS[ind, 2]
 207                 }
 208                 if (!is.null(tr$node.label)) {
 209                     for (j in 1:length(tr$node.label)) {
 210                         ind <- which(tr$node.label[j] == TRANS[, 1])
 211                         tr$node.label[j] <- TRANS[ind, 2]
 212                     }
 213                 }
 214                 trees[[i]] <- tr
 215             }
 216             translation <- FALSE
 217         }
 218     }
 219     for (i in 1:Ntree) {
 220         tr <- trees[[i]]
 221         ## Check here that the root edge is not incorrectly represented
 222         ## in the object of class "phylo" by simply checking that there
 223         ## is a bifurcation at the root
 224         if (!translation) n <- length(tr$tip.label)
 225         ROOT <- n + 1
 226         if (sum(tr$edge[, 1] == ROOT) == 1 && dim(tr$edge)[1] > 1) {
 227             stop(paste("The tree has apparently singleton node(s): cannot read tree file.\n  Reading NEXUS file aborted at tree no.", i, sep = ""))
 228         }
 229     }
 230     if (Ntree == 1) {
 231         trees <- trees[[1]]
 232         if (translation) {
 233             trees$tip.label <-
 234                 if (length(colon)) TRANS[, 2] else
 235                 TRANS[, 2][as.numeric(trees$tip.label)]
 236         }
 237     } else {
 238         if (!is.null(tree.names)) names(trees) <- tree.names
 239         if (translation) {
 240             if (length(colon) == Ntree) # .treeBuildWithTokens() was used
 241                 attr(trees, "TipLabel") <- TRANS[, 2]
 242             else { # reassign the tip labels then compress
 243                 for (i in 1:Ntree)
 244                     trees[[i]]$tip.label <-
 245                         TRANS[, 2][as.numeric(trees[[i]]$tip.label)]
 246                 trees <- .compressTipLabel(trees)
 247             }
 248         }
 249         class(trees) <- "multiPhylo"
 250         if (!all(nms.trees == "")) names(trees) <- nms.trees
 251     }
 252     trees
 253 }