R/rotate.R

   1 ### ROTATE
   2 ### Last update CH on 09.08.2007
   3
   4 # Contents:
   5 # 1. rotate
   6
   7 # 1: This function swops sister clades in a phylogenetic tree.
   8 # Branch lengths are considered.
   9 # Arguments:
  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.
  13 # Author: C.Heibl
  14
  15
  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]
  21                 while(max(x) > tips){
  22                         x <- x[x > tips]
  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]
  26                         }
  27                 x
  28                 }
  29 # function starts here
  30 # definitions
  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)
  44                 }
  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)
  48         }
  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]
  54         N <- N[N != node]
  55         if (length(N) > 2) N <- N[polytom]
  56         CLADE1 <- N[1]
  57         CLADE2 <- N[2]
  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)]
  65                         } else {
  66                         c1 <- G[,3][G[,2] == CLADE1]
  67                         c2 <- G[,3][G[,2] == CLADE1]
  68                         }
  69                 if (CLADE2 > nb.tips){
  70                         c3 <- G[,3][G[,2] == CLADE2]
  71                         c4 <- G[,3][G[,2] == max(CLADE22)]
  72                         } else {
  73                         c3 <- G[,3][G[,2] == CLADE2]
  74                         c4 <- G[,3][G[,2] == CLADE2]
  75                         }
  76
  77 # create new phy$edge and  phy$edge.length
  78 if (c2+1 == c3){
  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,])
  81                         if (with.br.length)
  82                         phy$edge.length <- c(phy$edge.length[c3:c4], phy$edge.length[c1:c2], phy$edge.length[(c4+1):nb.edges])
  83                 }
  84         if (c1 !=1 && c4 == nb.edges){
  85                 phy$edge <- rbind(phy$edge[1:(c1-1),], phy$edge[c3:c4,], phy$edge[c1:c2,])
  86                         if (with.br.length)
  87                         phy$edge.length <- c(phy$edge.length[1:(c1-1)], phy$edge.length[c3:c4], phy$edge.length[c1:c2])
  88                 }
  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,])
  91                         if (with.br.length)
  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])
  93                 }
  94         if (c1 ==1 && c4 == nb.edges){
  95                 phy$edge <- rbind(phy$edge[c3:c4,], phy$edge[c1:c2,])
  96                         if (with.br.length)
  97                         phy$edge.length <- c(phy$edge.length[c3:c4], phy$edge.length[c1:c2])
  98                 }
  99         }
 100 else {
 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,])
 103                         if (with.br.length)
 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])
 105                 }
 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,])
 108                         if (with.br.length)
 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])
 110                 }
 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,])
 113                         if (with.br.length)
 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])
 115                         }
 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,])
 118                         if (with.br.length)
 119                         phy$edge.length <- c(phy$edge.length[c3:c4], phy$edge.length[(c2+1):(c3-1)], phy$edge.length[c1:c2])
 120                 }
 121         }
 122         S <- write.tree(phy)
 123     phy <- if (!with.br.length) clado.build(S) else tree.build(S)
 124         phy
 125         }