YuLab-SMU / ggtree

:christmas_tree:Visualization and annotation of phylogenetic trees
https://yulab-smu.top/treedata-book/
833 stars 173 forks source link

angle not equal in equal-angle layout #513

Closed li-yq closed 2 years ago

li-yq commented 2 years ago

The equal-angle layout algorithm is not properly implemented.

https://github.com/YuLab-SMU/ggtree/blob/2aad9c49892caecd170518c189e0935d47c31364/R/tree-utilities.R#L95-L123

In line 107, total angle of current branch is calculated from end - start, but start is changed in line 122. So alpha will be incorrect after the first iteration, making branchs clustered in the plot.

Comparison of expected result and current implementation:

```R library(tidyverse) library(tidytree) library(ggtree) layoutEqualAngle <- function(model, branch.length = "branch.length"){ tree <- as.phylo(model) if (! is.null(tree$edge.length)) { if (anyNA(tree$edge.length)) { warning("'edge.length' contains NA values...\n## setting 'edge.length' to NULL automatically when plotting the tree...") tree$edge.length <- NULL } } if (is.null(tree$edge.length) || branch.length == "none") { tree <- set_branch_length_cladogram(tree) } N <- treeio::Nnode2(tree) brlen <- numeric(N) brlen[tree$edge[,2]] <- tree$edge.length root <- tidytree::rootnode(tree) ## Convert Phylo tree to data.frame. ## df <- as.data.frame.phylo_(tree) df <- as_tibble(model) %>% mutate(isTip = ! .data$node %in% .data$parent) ## NOTE: Angles (start, end, angle) are in half-rotation units (radians/pi or degrees/180) ## create and assign NA to the following fields. df$x <- 0 df$y <- 0 df$start <- 0 # Start angle of segment of subtree. df$end <- 0 # End angle of segment of subtree df$angle <- 0 # Orthogonal angle to beta for tip labels. ## Initialize root node position and angles. df[root, "x"] <- 0 df[root, "y"] <- 0 df[root, "start"] <- 0 # 0-degrees df[root, "end"] <- 2 # 360-degrees df[root, "angle"] <- 0 # Angle label. df$branch.length <- brlen[df$node] # for cladogram ## Get number of tips for each node in tree. ## nb.sp <- sapply(1:N, function(i) length(get.offspring.tip(tree, i))) ## self_include = TRUE to return itself if the input node is a tip nb.sp <- vapply(1:N, function(i) length(tidytree::offspring(tree, i, tiponly = TRUE, self_include = TRUE)), numeric(1)) ## Get list of node id's. nodes <- ggtree:::getNodes_by_postorder(tree) for(curNode in nodes) { ## Get number of tips for current node. curNtip <- nb.sp[curNode] ## Get array of child node indexes of current node. ## children <- getChild(tree, curNode) children <- treeio::child(tree, curNode) ## Get "start" and "end" angles of a segment for current node in the data.frame. start <- df[curNode, "start"] end <- df[curNode, "end"] cur_x = df[curNode, "x"] cur_y = df[curNode, "y"] total_angle = end - start for (child in children) { ## Get the number of tips for child node. ntip.child <- nb.sp[child] ## Calculated in half radians. ## alpha: angle of segment for i-th child with ntips_ij tips. ## alpha = (left_angle - right_angle) * (ntips_ij)/(ntips_current) alpha <- total_angle * ntip.child / curNtip ## beta = angle of line from parent node to i-th child. beta <- start + alpha / 2 length.child <- df[child, "branch.length"] ## update geometry of data.frame. ## Calculate (x,y) position of the i-th child node from current node. df[child, "x"] <- cur_x + cospi(beta) * length.child df[child, "y"] <- cur_y + sinpi(beta) * length.child ## Calculate orthogonal angle to beta for tip label. df[child, "angle"] <- -90 - 180 * beta * sign(beta - 1) ## Update the start and end angles of the childs segment. df[child, "start"] <- start df[child, "end"] <- start + alpha start <- start + alpha } } tree_df <- as_tibble(df) class(tree_df) <- c("tbl_tree", class(tree_df)) return(tree_df) } set.seed(20) tree <- rtree(20) tree |> layoutEqualAngle() |> ggplot() + geom_tree(layout = "slanted") + coord_fixed() tree |> ggtree:::layoutEqualAngle() |> ggplot() + geom_tree(layout = "slanted") + coord_fixed() ```

image

GuangchuangYu commented 2 years ago

thanks for your feedback. It was fixed in ggtree v >= 3.4.1.