2 ### Last update CH on 09.08.2007
7 # 1: This function swops sister clades in a phylogenetic tree.
8 # Branch lengths are considered.
10 # phy: an object of class phylo
11 # node: the number (integer) of the corresponding node or number or names of two tips that coalesce to th internal node
12 # polytom: use a vector of two integers to define those two clades of a tritomy, that are swopped. The default is c(1,2). The clade number is counted from the from bottom to top in the plotted tree.
16 rotate <- function(phy, node, polytom = c(1,2)){
17 # load DESCENDANTS function
18 DESCENDANTS <- function(tree, node){
19 tips <- length(tree$tip.label)
20 x <- tree$edge[,2][tree$edge[,1] == node]
23 for(h in 1:length(x)) tree$edge <- tree$edge[!tree$edge[,2] == x[h],]
24 for(i in 1:length(x)) tree$edge[,1][tree$edge[,1] == x[i]] <- node
25 x <- tree$edge[,2][tree$edge[,1] == node]
29 # function starts here
31 if (!inherits(phy, "phylo")) # is phy of class phylo?
32 stop("object \"phy\" is not of class \"phylo\"")
33 nb.tips <- length(phy$tip.label) # number of tiplabels
34 max.int.node <- phy$Nnode+nb.tips # number of last internal node
35 nb.edges <- dim(phy$edge)[1] # number of branches
36 if (length(node) == 2){ # get MRCA if tips are given for node
37 if (mode(node) == "character"){
38 if (any(!node %in% phy$tip.label)) # do tiplabels correspond
39 stop("object \"node\" contains tiplabels not present in object \"phy\"")
40 tips <- cbind(phy$tip.label, 1:nb.tips)
41 node[1] <- tips[,2][tips[,1] == node[1]]
42 node[2] <- tips[,2][tips[,1] == node[2]]
43 node <- as.numeric(node)
45 if (any(!node %in% 1:nb.tips)) # is phy of class phylo?
46 stop("object \"node\" does not contain terminal nodes")
47 node <- getMRCA(phy, node)
49 if (node <= nb.tips || node > max.int.node) # is node really internal?
50 stop("object \"node\" is not an internal node of object \"phy\"")
51 with.br.length <- !is.null(phy$edge.length) # does phy contain brlength?
52 G <- cbind(phy$edge, 1:(length(phy$edge)/2))
53 N <- phy$edge[phy$edge[,1] == node]
55 if (length(N) > 2) N <- N[polytom]
58 # do clades comprise interior nodes?
59 if (CLADE1 > nb.tips) CLADE11 <- DESCENDANTS(phy, CLADE1)
60 if (CLADE2 > nb.tips) CLADE22 <- DESCENDANTS(phy, CLADE2)
61 # calculate inidices of clades in phy.edge
62 if (CLADE1 > nb.tips){
63 c1 <- G[,3][G[,2] == CLADE1]
64 c2 <- G[,3][G[,2] == max(CLADE11)]
66 c1 <- G[,3][G[,2] == CLADE1]
67 c2 <- G[,3][G[,2] == CLADE1]
69 if (CLADE2 > nb.tips){
70 c3 <- G[,3][G[,2] == CLADE2]
71 c4 <- G[,3][G[,2] == max(CLADE22)]
73 c3 <- G[,3][G[,2] == CLADE2]
74 c4 <- G[,3][G[,2] == CLADE2]
77 # create new phy$edge and phy$edge.length
79 if (c1 == 1 && c4 != nb.edges){
80 phy$edge <- rbind(phy$edge[c3:c4,], phy$edge[c1:c2,], phy$edge[(c4+1):nb.edges,])
82 phy$edge.length <- c(phy$edge.length[c3:c4], phy$edge.length[c1:c2], phy$edge.length[(c4+1):nb.edges])
84 if (c1 !=1 && c4 == nb.edges){
85 phy$edge <- rbind(phy$edge[1:(c1-1),], phy$edge[c3:c4,], phy$edge[c1:c2,])
87 phy$edge.length <- c(phy$edge.length[1:(c1-1)], phy$edge.length[c3:c4], phy$edge.length[c1:c2])
89 if (c1 !=1 && c4 != nb.edges){
90 phy$edge <- rbind(phy$edge[1:(c1-1),], phy$edge[c3:c4,], phy$edge[c1:c2,], phy$edge[(c4+1):nb.edges,])
92 phy$edge.length <- c(phy$edge.length[1:(c1-1)], phy$edge.length[c3:c4], phy$edge.length[c1:c2], phy$edge.length[(c4+1):nb.edges])
94 if (c1 ==1 && c4 == nb.edges){
95 phy$edge <- rbind(phy$edge[c3:c4,], phy$edge[c1:c2,])
97 phy$edge.length <- c(phy$edge.length[c3:c4], phy$edge.length[c1:c2])
101 if (c1 == 1 && c4 != nb.edges){
102 phy$edge <- rbind(phy$edge[c3:c4,], phy$edge[(c2+1):(c3-1),], phy$edge[c1:c2,], phy$edge[(c4+1):nb.edges,])
104 phy$edge.length <- c(phy$edge.length[c3:c4], phy$edge.length[(c2+1):(c3-1)], phy$edge.length[c1:c2], phy$edge.length[(c4+1):nb.edges])
106 if (c1 !=1 && c4 == nb.edges){
107 phy$edge <- rbind(phy$edge[1:(c1-1),], phy$edge[c3:c4,], phy$edge[(c2+1):(c3-1),], phy$edge[c1:c2,])
109 phy$edge.length <- c(phy$edge.length[1:(c1-1)], phy$edge.length[c3:c4], phy$edge.length[(c2+1):(c3-1)], phy$edge.length[c1:c2])
111 if (c1 !=1 && c4 != nb.edges){
112 phy$edge <- rbind(phy$edge[1:(c1-1),], phy$edge[c3:c4,], phy$edge[(c2+1):(c3-1),], phy$edge[c1:c2,], phy$edge[(c4+1):nb.edges,])
114 phy$edge.length <- c(phy$edge.length[1:(c1-1)], phy$edge.length[c3:c4], phy$edge.length[(c2+1):(c3-1)], phy$edge.length[c1:c2], phy$edge.length[(c4+1):nb.edges])
116 if (c1 ==1 && c4 == nb.edges){
117 phy$edge <- rbind(phy$edge[c3:c4,], phy$edge[(c2+1):(c3-1),], phy$edge[c1:c2,])
119 phy$edge.length <- c(phy$edge.length[c3:c4], phy$edge.length[(c2+1):(c3-1)], phy$edge.length[c1:c2])
123 phy <- if (!with.br.length) clado.build(S) else tree.build(S)