R/read.nexus.R

   1 ## read.nexus.R (2007-12-22)
   2
   3 ##   Read Tree File in Nexus Format
   4
   5 ## Copyright 2003-2007 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 clado.build <- function(tp) {
  11     add.internal <- function() {
  12         edge[j, 1] <<- current.node
  13         node <<- node + 1
  14         edge[j, 2] <<- current.node <<- node
  15         j <<- j + 1
  16     }
  17     add.terminal <- function() {
  18         edge[j, 1] <<- current.node
  19         edge[j, 2] <<- tip
  20         tip.label[tip] <<- tpc[k]
  21         k <<- k + 1
  22         tip <<- tip + 1
  23         j <<- j + 1
  24     }
  25     go.down <- function() {
  26         l <- which(edge[, 2] == current.node)
  27         node.label[current.node - nb.tip] <<- tpc[k]
  28         k <<- k + 1
  29         current.node <<- edge[l, 1]
  30     }
  31     if (!length(grep(",", tp))) {
  32         obj <- list(edge = matrix(c(2, 1), 1, 2), Nnode = 1)
  33         tp <- unlist(strsplit(tp, "[\\(\\);]"))
  34         obj$tip.label <- tp[2]
  35         if (length(tp) == 3) obj$node.label <- tp[3]
  36         class(obj) <- "phylo"
  37         return(obj)
  38     }
  39     tsp <- unlist(strsplit(tp, NULL))
  40     tp <- gsub(")", ")NA", tp)
  41     tp <- gsub(" ", "", tp)
  42     tpc <- unlist(strsplit(tp, "[\\(\\),;]"))
  43     tpc <- tpc[tpc != ""]
  44     skeleton <- tsp[tsp == "(" | tsp == ")" | tsp == "," | tsp == ";"]
  45     nsk <- length(skeleton)
  46     nb.node <- length(skeleton[skeleton == ")"])
  47     nb.tip <- length(skeleton[skeleton == ","]) + 1
  48     ## We will assume there is an edge at the root;
  49     ## if so, it will be removed and put in a vector
  50     nb.edge <- nb.node + nb.tip
  51     node.label <- character(nb.node)
  52     tip.label <- character(nb.tip)
  53
  54     edge <- matrix(NA, nb.edge, 2)
  55     current.node <- node <- nb.tip + 1 # node number
  56     edge[nb.edge, 1] <- 0    # see comment above
  57     edge[nb.edge, 2] <- node #
  58
  59     ## j: index of the line number of edge
  60     ## k: index of the line number of tpc
  61     ## tip: tip number
  62     j <- k <- tip <- 1
  63
  64     for (i in 2:nsk) {
  65         if (skeleton[i] == "(") add.internal()      # add an internal branch (on top)
  66         if (skeleton[i] == ",") {
  67             if (skeleton[i - 1] != ")") add.terminal()   # add a terminal branch
  68         }
  69         if (skeleton[i] == ")") {
  70             if (skeleton[i - 1] == ",") {   # add a terminal branch and go down one level
  71                 add.terminal()
  72                 go.down()
  73             }
  74             if (skeleton[i - 1] == ")") go.down()   # go down one level
  75         }
  76     }
  77 #    if(node.label[1] == "NA") node.label[1] <- ""
  78     edge <- edge[-nb.edge, ]
  79     obj <- list(edge = edge, tip.label = tip.label,
  80                 Nnode = nb.node, node.label = node.label)
  81     obj$node.label <- if (all(obj$node.label == "NA")) NULL else gsub("^NA", "", obj$node.label)
  82     class(obj) <- "phylo"
  83     return(obj)
  84 }
  85
  86 read.nexus <- function(file, tree.names = NULL)
  87 {
  88     X <- scan(file = file, what = character(), sep = "\n", quiet = TRUE)
  89     ## first remove all the comments
  90     LEFT <- grep("\\[", X)
  91     RIGHT <- grep("\\]", X)
  92     if (length(LEFT)) {
  93         for (i in length(LEFT):1) {
  94             if (LEFT[i] == RIGHT[i]) {
  95                 X[LEFT[i]] <- gsub("\\[.*\\]", "", X[LEFT[i]])
  96             } else {
  97                 X[LEFT[i]] <- gsub("\\[.*", "", X[LEFT[i]])
  98                 X[RIGHT[i]] <- gsub(".*\\]", "", X[RIGHT[i]])
  99                 if (LEFT[i] < RIGHT[i] - 1) X <- X[-((LEFT[i] + 1):(RIGHT[i] - 1))]
 100             }
 101         }
 102     }
 103     X <- gsub("ENDBLOCK;", "END;", X, ignore.case = TRUE)
 104     endblock <- grep("END;", X, ignore.case = TRUE)
 105     semico <- grep(";", X)
 106     i1 <- grep("BEGIN TREES;", X, ignore.case = TRUE)
 107     i2 <- grep("TRANSLATE", X, ignore.case = TRUE)
 108     translation <- FALSE
 109     if (length(i2) == 1) if (i2 > i1) translation <- TRUE
 110     if (translation) {
 111         end <- semico[semico > i2][1]
 112         x <- paste(X[i2:end], sep = "", collapse = "")
 113         x <- gsub("TRANSLATE", "", x, ignore.case = TRUE)
 114         x <- unlist(strsplit(x, "[,; \t]"))
 115         x <- x[x != ""]
 116         TRANS <- matrix(x, ncol = 2, byrow = TRUE)
 117         TRANS[, 2] <- gsub("['\"]", "", TRANS[, 2])
 118     }
 119     start <- if (translation)  semico[semico > i2][1] + 1 else semico[semico > i1][1]
 120     end <- endblock[endblock > i1][1] - 1
 121     tree <- paste(X[start:end], sep = "", collapse = "")
 122     tree <- gsub(" ", "", tree)
 123     tree <- unlist(strsplit(tree, "[=;]"))
 124     tree <- tree[grep("[\\(\\)]", tree)]
 125     nb.tree <- length(tree)
 126     STRING <- as.list(tree)
 127     trees <- list()
 128     for (i in 1:nb.tree) {
 129         obj <- if (length(grep(":", STRING[[i]]))) tree.build(STRING[[i]]) else clado.build(STRING[[i]])
 130         if (translation) {
 131             for (j in 1:length(obj$tip.label)) {
 132                 ind <- which(obj$tip.label[j] == TRANS[, 1])
 133                 obj$tip.label[j] <- TRANS[ind, 2]
 134             }
 135             if (!is.null(obj$node.label)) {
 136                 for (j in 1:length(obj$node.label)) {
 137                     ind <- which(obj$node.label[j] == TRANS[, 1])
 138                     obj$node.label[j] <- TRANS[ind, 2]
 139                 }
 140             }
 141         }
 142         ## Check here that the root edge is not incorrectly represented
 143         ## in the object of class "phylo" by simply checking that there
 144         ## is a bifurcation at the root
 145         ROOT <- length(obj$tip.label) + 1
 146         if (sum(obj$edge[, 1] == ROOT) == 1 && dim(obj$edge)[1] > 1) {
 147             stop(paste("There is apparently two root edges in your file: cannot read tree file.\n  Reading NEXUS file aborted at tree no.", i, sep = ""))
 148         }
 149         trees[[i]] <- obj
 150     }
 151     if (nb.tree == 1) trees <- trees[[1]] else {
 152         names(trees) <- if (is.null(tree.names))
 153             paste("tree", 1:nb.tree, sep = "") else tree.names
 154         class(trees) <- "multiPhylo"
 155     }
 156     if (length(grep("[\\/]", file)) == 1) attr(trees, "origin") <- file
 157     else attr(trees, "origin") <- paste(getwd(), file, sep = "/")
 158     trees
 159 }