1 ## read.nexus.R (2007-12-22)
3 ## Read Tree File in Nexus Format
5 ## Copyright 2003-2007 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 clado.build <- function(tp) {
11 add.internal <- function() {
12 edge[j, 1] <<- current.node
14 edge[j, 2] <<- current.node <<- node
17 add.terminal <- function() {
18 edge[j, 1] <<- current.node
20 tip.label[tip] <<- tpc[k]
25 go.down <- function() {
26 l <- which(edge[, 2] == current.node)
27 node.label[current.node - nb.tip] <<- tpc[k]
29 current.node <<- edge[l, 1]
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]
39 tsp <- unlist(strsplit(tp, NULL))
40 tp <- gsub(")", ")NA", tp)
41 tp <- gsub(" ", "", tp)
42 tpc <- unlist(strsplit(tp, "[\\(\\),;]"))
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)
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 #
59 ## j: index of the line number of edge
60 ## k: index of the line number of tpc
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
69 if (skeleton[i] == ")") {
70 if (skeleton[i - 1] == ",") { # add a terminal branch and go down one level
74 if (skeleton[i - 1] == ")") go.down() # go down one level
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)
86 read.nexus <- function(file, tree.names = NULL)
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)
93 for (i in length(LEFT):1) {
94 if (LEFT[i] == RIGHT[i]) {
95 X[LEFT[i]] <- gsub("\\[.*\\]", "", X[LEFT[i]])
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))]
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)
109 if (length(i2) == 1) if (i2 > i1) translation <- TRUE
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]"))
116 TRANS <- matrix(x, ncol = 2, byrow = TRUE)
117 TRANS[, 2] <- gsub("['\"]", "", TRANS[, 2])
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)
128 for (i in 1:nb.tree) {
129 obj <- if (length(grep(":", STRING[[i]]))) tree.build(STRING[[i]]) else clado.build(STRING[[i]])
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]
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]
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 = ""))
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"
156 if (length(grep("[\\/]", file)) == 1) attr(trees, "origin") <- file
157 else attr(trees, "origin") <- paste(getwd(), file, sep = "/")