@@ -27,6 +27,7 @@ Imports:

     scales,

     tibble,

     tidyr,

+    tidytree,

     utils

 Suggests:

     Biostrings,

@@ -7,7 +7,8 @@ all: rd readme check clean

 alldocs: site rd readme

 rd:

-	Rscript -e 'roxygen2::roxygenise(".")'

+	Rscript -e 'library(methods); devtools::document()'

+# Rscript -e 'roxygen2::roxygenise(".")'

 readme:

 	Rscript -e 'rmarkdown::render("README.Rmd")'

@@ -28,11 +29,11 @@ install:

 	cd ..;\

 	R CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz

-check: build

+check: rd build

 	cd ..;\

 	Rscript -e 'rcmdcheck::rcmdcheck("$(PKGNAME)_$(PKGVERS).tar.gz")'

-check2: build

+check2: rd build

 	cd ..;\

 	R CMD check $(PKGNAME)_$(PKGVERS).tar.gz

@@ -2,8 +2,6 @@

 S3method(as.binary,phylo)

 S3method(as.data.frame,phylo)

-S3method(as_data_frame,phylo)

-S3method(as_data_frame,treedata)

 S3method(fortify,multiPhylo)

 S3method(fortify,obkData)

 S3method(fortify,phylo)

@@ -184,13 +182,14 @@ importFrom(methods,setGeneric)

 importFrom(methods,setOldClass)

 importFrom(rvcheck,get_fun_from_pkg)

 importFrom(scales,alpha)

-importFrom(tibble,as_data_frame)

 importFrom(tibble,data_frame)

 importFrom(tidyr,gather)

+importFrom(tidytree,as_data_frame)

 importFrom(treeio,Nnode)

 importFrom(treeio,Ntip)

 importFrom(treeio,as.phylo)

 importFrom(treeio,as.treedata)

+importFrom(treeio,get_tree_data)

 importFrom(treeio,groupClade)

 importFrom(treeio,groupOTU)

 importFrom(utils,modifyList)

new file mode 100644

@@ -0,0 +1,63 @@

+## not used.

+## we don't guess mrsd from tip labels

+## user should provide mrsd parameter in ggtree if they want to use time-scaled tree

+scaleX_by_time <- function(df, as.Date=FALSE) {

+    time <- with(df, gsub(".*[_/]{1}(\\d+\\.*\\d+)$", "\\1", label[isTip])) %>% as.numeric

+    latest <- which.max(time)

+

+    scaleX_by_time_from_mrsd(df, decimal2Date(time[latest]), as.Date)

+}

+

+

+scaleX_by_time_from_mrsd <- function(df, mrsd, as.Date) {

+    mrsd %<>% as.Date

+    date <- Date2decimal(mrsd)

+

+    df$x <- df$x + date - max(df$x)

+    df$branch <- with(df, (x[match(parent, node)] + x)/2)

+    ## df$branch <- (df[df$parent, "x"] + df[, "x"])/2

+

+    if (as.Date) {

+        df$x <- decimal2Date(df$x)

+        df$branch <- decimal2Date(df$branch)

+    }

+

+    return(df)

+}

+

+

+##' convert Date to decimal format, eg "2014-05-05" to "2014.34"

+##'

+##'

+##' @title Date2decimal

+##' @param x Date

+##' @return numeric

+##' @export

+##' @author Guangchuang Yu

+Date2decimal <- function(x) {

+    if (is(x, "numeric")) {

+        return(x)

+    }

+

+    if (is(x, "character")) {

+        x <- as.Date(x)

+    }

+    year <- format(x, "%Y")

+    y <- x - as.Date(paste0(year, "-01-01"))

+    as.numeric(year) + as.numeric(y)/365

+}

+

+##' convert decimal format to Date, eg "2014.34" to "2014-05-05"

+##'

+##'

+##' @title decimal2Date

+##' @param x numerical number, eg 2014.34

+##' @return Date

+##' @export

+##' @author Guangchuang Yu

+decimal2Date <- function(x) {

+    date <- as.Date(paste0(floor(x), "-01-01"))

+    date + as.numeric(sub("^\\d+", "0", x)) * 365

+}

+

+

@@ -8,6 +8,7 @@

 ##' @return polytomy tree

 ##' @author Guangchuang

 ##' @importFrom ape di2multi

+##' @importFrom treeio Ntip

 ##' @export

 as.polytomy <- function(tree, feature, fun) {

     if (!is(tree, 'phylo')) {

new file mode 100644

@@ -0,0 +1,158 @@

+##' @importFrom ggplot2 fortify

+##' @method fortify treedata

+##' @export

+fortify.treedata <- function(model, data,

+                             layout        = "rectangular",

+                             yscale        = "none",

+                             ladderize     = TRUE,

+                             right         = FALSE,

+                             branch.length = "branch.length",

+                             mrsd          = NULL,

+                             as.Date       = FALSE, ...) {

+

+    model <- set_branch_length(model, branch.length)

+

+    fortify.phylo(model, data,

+                  layout        = layout,

+                  yscale        = yscale,

+                  ladderize     = ladderize,

+                  right         = right,

+                  branch.length = branch.length,

+                  mrsd          = mrsd,

+                  as.Date       = as.Date, ...)

+}

+

+##' @importFrom ape ladderize

+##' @importFrom treeio as.phylo

+##' @importFrom treeio Nnode

+##' @importFrom tibble data_frame

+##' @importFrom dplyr full_join

+##' @importFrom tidytree as_data_frame

+##' @method fortify phylo

+##' @export

+fortify.phylo <- function(model, data,

+                          layout        = "rectangular",

+                          ladderize     = TRUE,

+                          right         = FALSE,

+                          branch.length = "branch.length",

+                          mrsd          = NULL,

+                          as.Date       = FALSE,

+                          yscale        = "none",

+                          ...) {

+

+    x <- as.phylo(model) ## reorder.phylo(get.tree(model), "postorder")

+    if (ladderize == TRUE) {

+        x <- ladderize(x, right=right)

+    }

+

+    if (! is.null(x$edge.length)) {

+        if (anyNA(x$edge.length)) {

+            warning("'edge.length' contains NA values...\n## setting 'edge.length' to NULL automatically when plotting the tree...")

+            x$edge.length <- NULL

+        }

+    }

+

+    if (is.null(x$edge.length) || branch.length == "none") {

+        xpos <- getXcoord_no_length(x)

+    } else {

+        xpos <- getXcoord(x)

+    }

+

+    ypos <- getYcoord(x)

+    N <- Nnode(x, internal.only=FALSE)

+    xypos <- data_frame(node=1:N, x=xpos, y=ypos)

+

+    df <- as_data_frame(model)

+

+    res <- full_join(df, xypos, by = "node")

+

+    ## add branch mid position

+    res <- calculate_branch_mid(res)

+

+    if (!is.null(mrsd)) {

+        res <- scaleX_by_time_from_mrsd(res, mrsd, as.Date)

+    }

+

+    if (layout == "slanted") {

+        res <- add_angle_slanted(res)

+    } else {

+        ## angle for all layout, if 'rectangular', user use coord_polar, can still use angle

+        res <- calculate_angle(res)

+    }

+    scaleY(as.phylo(model), res, yscale, layout, ...)

+}

+

+##' @importFrom treeio get_tree_data

+set_branch_length <- function(tree_object, branch.length) {

+    if (branch.length == "branch.length") {

+        return(tree_object)

+    } else if (branch.length == "none") {

+        tree_object@phylo$edge.length <- NULL

+        return(tree_object)

+    }

+

+    if (is(tree_object, "phylo")) {

+        return(tree_object)

+    }

+

+    tree_anno <- get_tree_data(tree_object)

+    tree_anno$node <- as.integer(tree_anno$node)

+

+    phylo <- as.phylo(tree_object)

+

+    cn <- colnames(tree_anno)

+    cn <- cn[!cn %in% c('node', 'parent')]

+

+    length <- match.arg(branch.length, cn)

+

+    if (all(is.na(as.numeric(tree_anno[[length]])))) {

+        stop("branch.length should be numerical attributes...")

+    }

+

+    edge <- as_data_frame(phylo$edge)

+    colnames(edge) <- c("parent", "node")

+

+    dd <- full_join(edge, tree_anno, by = "node")

+

+    dd <- dd[match(edge[['node']], dd[['node']]),]

+    len <- unlist(dd[[length]])

+    len <- as.numeric(len)

+    len[is.na(len)] <- 0

+

+    phylo$edge.length <- len

+

+    tree_object@phylo <- phylo

+    return(tree_object)

+}

+

+

+calculate_angle <- function(data) {

+    data$angle <- 360/(diff(range(data$y)) + 1) * data$y

+    return(data)

+}

+

+

+

+scaleY <- function(phylo, df, yscale, layout, ...) {

+    if (yscale == "none") {

+        return(df)

+    }

+    if (! yscale %in% colnames(df)) {

+        warning("yscale is not available...\n")

+        return(df)

+    }

+    if (is.numeric(df[[yscale]])) {

+        y <- getYcoord_scale_numeric(phylo, df, yscale, ...)

+        ## if (order.y) {

+        ##     y <- getYcoord_scale2(phylo, df, yscale)

+        ## } else {

+        ##     y <- getYcoord_scale(phylo, df, yscale)

+        ## }

+    } else {

+        y <- getYcoord_scale_category(phylo, df, yscale, ...)

+    }

+

+    df[, "y"] <- y

+

+    return(df)

+}

@@ -17,8 +17,6 @@

 inset <- function(tree_view, insets, width=0.1, height=0.1, hjust=0, vjust=0,

                   x="node", reverse_x=FALSE, reverse_y=FALSE) {

-    message("The inset function will be defunct in next release, please use ggimage::geom_subview() instead.")

-

     df <- tree_view$data[as.numeric(names(insets)),]

     x <- match.arg(x, c("node", "branch", "edge"))

@@ -764,3 +764,717 @@ getNodesBreadthFirst.df <- function(df){

+##' convert tip or node label(s) to internal node number

+##'

+##'

+##' @title nodeid

+##' @param x tree object or graphic object return by ggtree

+##' @param label tip or node label(s)

+##' @return internal node number

+##' @importFrom methods is

+##' @export

+##' @author Guangchuang Yu

+nodeid <- function(x, label) {

+    if (is(x, "gg"))

+        return(nodeid.gg(x, label))

+

+    nodeid.tree(x, label)

+}

+

+nodeid.tree <- function(tree, label) {

+    tr <- get.tree(tree)

+    lab <- c(tr$tip.label, tr$node.label)

+    match(label, lab)

+}

+

+nodeid.gg <- function(p, label) {

+    p$data$node[match(label, p$data$label)]

+}

+

+

+reroot_node_mapping <- function(tree, tree2) {

+    root <- getRoot(tree)

+

+    node_map <- data.frame(from=1:getNodeNum(tree), to=NA, visited=FALSE)

+    node_map[1:Ntip(tree), 2] <- match(tree$tip.label, tree2$tip.label)

+    node_map[1:Ntip(tree), 3] <- TRUE

+

+    node_map[root, 2] <- root

+    node_map[root, 3] <- TRUE

+

+    node <- rev(tree$edge[,2])

+    for (k in node) {

+        ip <- getParent(tree, k)

+        if (node_map[ip, "visited"])

+            next

+

+        cc <- getChild(tree, ip)

+        node2 <- node_map[cc,2]

+        if (anyNA(node2)) {

+            node <- c(node, k)

+            next

+        }

+

+        to <- unique(sapply(node2, getParent, tr=tree2))

+        to <- to[! to %in% node_map[,2]]

+        node_map[ip, 2] <- to

+        node_map[ip, 3] <- TRUE

+    }

+    node_map <- node_map[, -3]

+    return(node_map)

+}

+

+

+

+##' Get parent node id of child node.

+##'

+##' @title getParent.df

+##' @param df tree data.frame

+##' @param node is the node id of child in tree.

+##' @return integer node id of parent

+getParent.df <- function(df, node) {

+    i <- which(df$node == node)

+    parent_id <- df$parent[i]

+    if (parent_id == node | is.na(parent_id)) {

+        ## root node

+        return(0)

+    }

+    return(parent_id)

+}

+

+

+getAncestor.df <- function(df, node) {

+    anc <- getParent.df(df, node)

+    anc <- anc[anc != 0]

+    if (length(anc) == 0) {

+        # stop("selected node is root...")

+      return(0)

+    }

+    i <- 1

+    while(i<= length(anc)) {

+        anc <- c(anc, getParent.df(df, anc[i]))

+        anc <- anc[anc != 0]

+        i <- i+1

+    }

+    return(anc)

+}

+

+

+

+##' Get list of child node id numbers of parent node

+##'

+##' @title getChild.df

+##' @param df tree data.frame

+##' @param node is the node id of child in tree.

+##' @return list of child node ids of parent

+getChild.df <- function(df, node) {

+    i <- which(df$parent == node)

+    if (length(i) == 0) {

+        return(0) # it has no children, hence tip node.

+    }

+    res <- df$node[i]

+    res <- res[res != node] ## node may root

+    return(res)

+}

+

+get.offspring.df <- function(df, node) {

+    sp <- getChild.df(df, node)

+    sp <- sp[sp != 0] # Remove root node.

+    if (length(sp) == 0) {

+        #stop("input node is a tip...")

+      return(0)

+    }

+

+    i <- 1

+    while(i <= length(sp)) {

+        sp <- c(sp, getChild.df(df, sp[i]))

+        sp <- sp[sp != 0]

+        i <- i + 1

+    }

+    return(sp)

+}

+

+

+

+##' extract offspring tips

+##'

+##'

+##' @title get.offspring.tip

+##' @param tr tree

+##' @param node node

+##' @return tip label

+##' @author ygc

+##' @importFrom ape extract.clade

+##' @export

+get.offspring.tip <- function(tr, node) {

+    if ( ! node %in% tr$edge[,1]) {

+        ## return itself

+        return(tr$tip.label[node])

+    }

+    clade <- extract.clade(tr, node)

+    clade$tip.label

+}

+

+

+

+

+getParent <- function(tr, node) {

+    if ( node == getRoot(tr) )

+        return(0)

+    edge <- tr[["edge"]]

+    parent <- edge[,1]

+    child <- edge[,2]

+    res <- parent[child == node]

+    if (length(res) == 0) {

+        stop("cannot found parent node...")

+    }

+    if (length(res) > 1) {

+        stop("multiple parent found...")

+    }

+    return(res)

+}

+

+

+

+

+getChild <- function(tr, node) {

+    # Get edge matrix from phylo object.

+    edge <- tr[["edge"]]

+    # Select all rows that match "node".

+    res <- edge[edge[,1] == node, 2]

+    ## if (length(res) == 0) {

+    ##     ## is a tip

+    ##     return(NA)

+    ## }

+    return(res)

+}

+

+

+getSibling <- function(tr, node) {

+    root <- getRoot(tr)

+    if (node == root) {

+        return(NA)

+    }

+

+    parent <- getParent(tr, node)

+    child <- getChild(tr, parent)

+    sib <- child[child != node]

+    return(sib)

+}

+

+

+getAncestor <- function(tr, node) {

+    root <- getRoot(tr)

+    if (node == root) {

+        return(NA)

+    }

+    parent <- getParent(tr, node)

+    res <- parent

+    while(parent != root) {

+        parent <- getParent(tr, parent)

+        res <- c(res, parent)

+    }

+    return(res)

+}

+

+

+isRoot <- function(tr, node) {

+    getRoot(tr) == node

+}

+

+getNodeName <- function(tr) {

+    if (is.null(tr$node.label)) {

+        n <- length(tr$tip.label)

+        nl <- (n + 1):(2 * n - 2)

+        nl <- as.character(nl)

+    }

+    else {

+        nl <- tr$node.label

+    }

+    nodeName <- c(tr$tip.label, nl)

+    return(nodeName)

+}

+

+

+

+get.trunk <- function(tr) {

+    root <- getRoot(tr)

+    path_length <- sapply(1:(root-1), function(x) get.path_length(tr, root, x))

+    i <- which.max(path_length)

+    return(get.path(tr, root, i))

+}

+

+##' path from start node to end node

+##'

+##'

+##' @title get.path

+##' @param phylo phylo object

+##' @param from start node

+##' @param to end node

+##' @return node vectot

+##' @export

+##' @author Guangchuang Yu

+get.path <- function(phylo, from, to) {

+    anc_from <- getAncestor(phylo, from)

+    anc_from <- c(from, anc_from)

+    anc_to <- getAncestor(phylo, to)

+    anc_to <- c(to, anc_to)

+    mrca <- intersect(anc_from, anc_to)[1]

+

+    i <- which(anc_from == mrca)

+    j <- which(anc_to == mrca)

+

+    path <- c(anc_from[1:i], rev(anc_to[1:(j-1)]))

+    return(path)

+}

+

+

+get.path_length <- function(phylo, from, to, weight=NULL) {

+    path <- get.path(phylo, from, to)

+    if (is.null(weight)) {

+        return(length(path)-1)

+    }

+

+    df <- fortify(phylo)

+    if ( ! (weight %in% colnames(df))) {

+        stop("weight should be one of numerical attributes of the tree...")

+    }

+

+    res <- 0

+

+    get_edge_index <- function(df, from, to) {

+        which((df[,1] == from | df[,2] == from) &

+                  (df[,1] == to | df[,2] == to))

+    }

+

+    for(i in 1:(length(path)-1)) {

+        ee <- get_edge_index(df, path[i], path[i+1])

+        res <- res + df[ee, weight]

+    }

+

+    return(res)

+}

+

+##' @importFrom ape reorder.phylo

+getNodes_by_postorder <- function(tree) {

+  tree <- reorder.phylo(tree, "postorder")

+    unique(rev(as.vector(t(tree$edge[,c(2,1)]))))

+}

+

+getXcoord2 <- function(x, root, parent, child, len, start=0, rev=FALSE) {

+    x[root] <- start

+    x[-root] <- NA  ## only root is set to start, by default 0

+

+    currentNode <- root

+    direction <- 1

+    if (rev == TRUE) {

+        direction <- -1

+    }

+

+    while(anyNA(x)) {

+        idx <- which(parent %in% currentNode)

+        newNode <- child[idx]

+        x[newNode] <- x[parent[idx]]+len[idx] * direction

+        currentNode <- newNode

+    }

+

+    return(x)

+}

+

+

+

+

+

+

+getXcoord_no_length <- function(tr) {

+    edge <- tr$edge

+    parent <- edge[,1]

+    child <- edge[,2]

+    root <- getRoot(tr)

+

+    len <- tr$edge.length

+

+    N <- getNodeNum(tr)

+    x <- numeric(N)

+    ntip <- Ntip(tr)

+    currentNode <- 1:ntip

+    x[-currentNode] <- NA

+

+    cl <- split(child, parent)

+    child_list <- list()

+    child_list[as.numeric(names(cl))] <- cl

+

+    while(anyNA(x)) {

+        idx <- match(currentNode, child)

+        pNode <- parent[idx]

+        ## child number table

+        p1 <- table(parent[parent %in% pNode])

+        p2 <- table(pNode)

+        np <- names(p2)

+        i <- p1[np] == p2

+        newNode <- as.numeric(np[i])

+

+        exclude <- rep(NA, max(child))

+        for (j in newNode) {

+            x[j] <- min(x[child_list[[j]]]) - 1

+            exclude[child_list[[j]]] <- child_list[[j]]

+        }

+        exclude <- exclude[!is.na(exclude)]

+

+        ## currentNode %<>% `[`(!(. %in% exclude))

+        ## currentNode %<>% c(., newNode) %>% unique

+        currentNode <- currentNode[!currentNode %in% exclude]

+        currentNode <- unique(c(currentNode, newNode))

+

+    }

+    x <- x - min(x)

+    return(x)

+}

+

+

+

+

+getXcoord <- function(tr) {

+    edge <- tr$edge

+    parent <- edge[,1]

+    child <- edge[,2]

+    root <- getRoot(tr)

+

+    len <- tr$edge.length

+

+    N <- getNodeNum(tr)

+    x <- numeric(N)

+    x <- getXcoord2(x, root, parent, child, len)

+    return(x)

+}

+

+

+

+

+## scale the branch (the line plotted) to the actual value of edge length

+## but it seems not the good idea as if we want to add x-axis (e.g. time-scaled tree)

+## then the x-value is not corresponding to edge length as in rectangular layout

+## getXYcoord_slanted <- function(tr) {

+##     edge <- tr$edge

+##     parent <- edge[,1]

+##     child <- edge[,2]

+##     root <- getRoot(tr)

+

+##     N <- getNodeNum(tr)

+##     len <- tr$edge.length

+##     y <- getYcoord(tr, step=min(len)/2)

+##     len <- sqrt(len^2 - (y[parent]-y[child])^2)

+##     x <- numeric(N)

+##     x <- getXcoord2(x, root, parent, child, len)

+##     res <- data.frame(x=x, y=y)

+##     return(res)

+## }

+

+

+

+## @importFrom magrittr %>%

+##' @importFrom magrittr equals

+getYcoord <- function(tr, step=1) {

+    Ntip <- length(tr[["tip.label"]])

+    N <- getNodeNum(tr)

+

+    edge <- tr[["edge"]]

+    parent <- edge[,1]

+    child <- edge[,2]

+

+    cl <- split(child, parent)

+    child_list <- list()

+    child_list[as.numeric(names(cl))] <- cl

+

+    y <- numeric(N)

+    tip.idx <- child[child <= Ntip]

+    y[tip.idx] <- 1:Ntip * step

+    y[-tip.idx] <- NA

+

+    ## use lookup table

+    pvec <- integer(max(tr$edge))

+    pvec[child] = parent

+

+    currentNode <- 1:Ntip

+    while(anyNA(y)) {

+        ## pNode <- unique(parent[child %in% currentNode])

+        pNode <- unique(pvec[currentNode])

+

+        ## piping of magrittr is slower than nested function call.

+        ## pipeR is fastest, may consider to use pipeR

+        ##

+        ## child %in% currentNode %>% which %>% parent[.] %>% unique

+        ## idx <- sapply(pNode, function(i) all(child[parent == i] %in% currentNode))

+        idx <- sapply(pNode, function(i) all(child_list[[i]] %in% currentNode))

+        newNode <- pNode[idx]

+

+        y[newNode] <- sapply(newNode, function(i) {

+            mean(y[child_list[[i]]], na.rm=TRUE)

+            ##child[parent == i] %>% y[.] %>% mean(na.rm=TRUE)

+        })

+

+        currentNode <- c(currentNode[!currentNode %in% unlist(child_list[newNode])], newNode)

+        ## currentNode <- c(currentNode[!currentNode %in% child[parent %in% newNode]], newNode)

+        ## parent %in% newNode %>% child[.] %>%

+        ##     `%in%`(currentNode, .) %>% `!` %>%

+        ##         currentNode[.] %>% c(., newNode)

+    }

+

+    return(y)

+}

+

+

+getYcoord_scale <- function(tr, df, yscale) {

+

+    N <- getNodeNum(tr)

+    y <- numeric(N)

+

+    root <- getRoot(tr)

+    y[root] <- 0

+    y[-root] <- NA

+

+    edge <- tr$edge

+    parent <- edge[,1]

+    child <- edge[,2]

+

+    currentNodes <- root

+    while(anyNA(y)) {

+        newNodes <- c()

+        for (currentNode in currentNodes) {

+            idx <- which(parent %in% currentNode)

+            newNode <- child[idx]

+            direction <- -1

+            for (i in seq_along(newNode)) {

+                y[newNode[i]] <- y[currentNode] + df[newNode[i], yscale] * direction

+                direction <- -1 * direction

+            }

+            newNodes <- c(newNodes, newNode)

+        }

+        currentNodes <- unique(newNodes)

+    }

+    if (min(y) < 0) {

+        y <- y + abs(min(y))

+    }

+    return(y)

+}

+

+

+getYcoord_scale2 <- function(tr, df, yscale) {

+    root <- getRoot(tr)

+

+    pathLength <- sapply(1:length(tr$tip.label), function(i) {

+        get.path_length(tr, i, root, yscale)

+    })

+

+    ordered_tip <- order(pathLength, decreasing = TRUE)

+    ii <- 1

+    ntip <- length(ordered_tip)

+    while(ii < ntip) {

+        sib <- getSibling(tr, ordered_tip[ii])

+        if (length(sib) == 0) {

+            ii <- ii + 1

+            next

+        }

+        jj <- which(ordered_tip %in% sib)

+        if (length(jj) == 0) {

+            ii <- ii + 1

+            next

+        }

+        sib <- ordered_tip[jj]

+        ordered_tip <- ordered_tip[-jj]

+        nn <- length(sib)

+        if (ii < length(ordered_tip)) {

+            ordered_tip <- c(ordered_tip[1:ii],sib, ordered_tip[(ii+1):length(ordered_tip)])

+        } else {

+            ordered_tip <- c(ordered_tip[1:ii],sib)

+        }

+

+        ii <- ii + nn + 1

+    }

+

+

+    long_branch <- getAncestor(tr, ordered_tip[1]) %>% rev

+    long_branch <- c(long_branch, ordered_tip[1])

+

+    N <- getNodeNum(tr)

+    y <- numeric(N)

+

+    y[root] <- 0

+    y[-root] <- NA

+

+    ## yy <- df[, yscale]

+    ## yy[is.na(yy)] <- 0

+

+    for (i in 2:length(long_branch)) {

+        y[long_branch[i]] <- y[long_branch[i-1]] + df[long_branch[i], yscale]

+    }

+

+    parent <- df[, "parent"]

+    child <- df[, "node"]

+

+    currentNodes <- root

+    while(anyNA(y)) {

+        newNodes <- c()

+        for (currentNode in currentNodes) {

+            idx <- which(parent %in% currentNode)

+            newNode <- child[idx]