@@ -1,11 +1,14 @@

 ##' drawing phylogenetic trees from list of phylo objects

+##' 

+##' The trees plotted by 'ggdensitree()' will be stacked on top of each other and the 

+##' structures of the trees will be rotated to ensure the consistency of the tip order. 

 ##'

 ##' @title ggdensitree

 ##' @param data a list of phylo objects or any object with an as.phylo and fortify method

 ##' @param mapping aesthetic mapping

 ##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')

 ##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; 'mode' to order the tips by the most common order; 'mds' to order the tips based on MDS of the path length between the tips;  or 'mds_dist' to order the tips based on MDS of the distance between the tips (default: 'mode')

-##' @param align.tips TRUE to align trees by their tips and FALSE to align trees by their root (default: TRUE)

+##' @param align.tips TRUE (default) to align trees by their tips and FALSE to align trees by their root

 ##' @param jitter deviation to jitter tips

 ##' @param ... additional parameters passed to fortify, ggtree and geom_tree

 ##' @return tree layer

@@ -13,6 +16,10 @@

 ##' @importFrom magrittr add

 ##' @export

 ##' @author Yu Guangchuang, Bradley R. Jones

+##' @references

+##' For more detailed demonstration of this function, please refer to chapter 4.4.2 of 

+##' *Data Integration, Manipulation and Visualization of Phylogenetic Trees*

+##' <http://yulab-smu.top/treedata-book/index.html> by Guangchuang Yu.

 ##' @examples

 ##' require(ape)

 ##' require(dplyr)

@@ -119,11 +119,10 @@ ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mo

 	## reorder tips (and shift x id align tips)

 	max.x <- vapply(trees.f, function(x) max(x$x, na.rm = TRUE), numeric(1))

-	farthest <- max(max.x)

 	trees.f <- lapply(1:length(trees), function(i) {

 		trees.f[[i]]$y <- getYcoord(trees[[i]], tip.order = tip.order)

 		if (align.tips) {

-			trees.f[[i]]$x <- trees.f[[i]]$x + (farthest - max.x[i])

+			trees.f[[i]]$x <- trees.f[[i]]$x - max.x[i]

 		}

 		if (i > 1 && jitter > 0) {

 			trees.f[[i]]$y[1:n.tips] %<>% add(stats::rnorm(n.tips, mean=0, sd=jitter))

@@ -62,15 +62,15 @@ ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mo

 	if (length(tip.order) == 1) {

 	  if (tip.order == 'mode') {

 	    first.label <- trees.f[[1]] %>%

-	      dplyr::filter(isTip) %>%

-	      dplyr::arrange(y) %>%

-	      dplyr::pull(label)

+	      dplyr::filter(.data$isTip) %>%

+	      dplyr::arrange(.data$y) %>%

+	      dplyr::pull(.data$label)

 		  res <- lapply(

 	      trees.f,

 	      . %>%

-	        dplyr::filter(isTip) %>%

-	        dplyr::arrange(y) %>%

+	        dplyr::filter(.data$isTip) %>%

+	        dplyr::arrange(.data$y) %>%

 	        dplyr::pull("label") %>%

 	        match(first.label)

 	    ) %>%

@@ -78,8 +78,8 @@ ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mo

 	      as.data.frame() %>%

 	      dplyr::group_by_all() %>%

 	      dplyr::summarize(count = dplyr::n(), .groups = 'drop') %>%

-	      dplyr::filter(count == max(count)) %>%

-	      dplyr::select(-count) %>%

+	      dplyr::filter(.data$count == max(.data$count)) %>%

+	      dplyr::select(-.data$count) %>%

 	      dplyr::slice(1) %>%

 	      unlist()

@@ -144,4 +144,4 @@ cophenetic.phylo.check.length <- function(tree, method) {

 		tree$edge.length <- rep(1, nrow(tree$edge))

 	stats::cophenetic(tree)

-}

\ No newline at end of file

+}

@@ -71,7 +71,7 @@ ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mo

 	      . %>%

 	        dplyr::filter(isTip) %>%

 	        dplyr::arrange(y) %>%

-	        `$`("label") %>%

+	        dplyr::pull("label") %>%

 	        match(first.label)

 	    ) %>%

 	      do.call(rbind, .) %>%

@@ -4,7 +4,7 @@

 ##' @param data a list of phylo objects or any object with an as.phylo and fortify method

 ##' @param mapping aesthetic mapping

 ##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')

-##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; 'mds' to order the tips based on MDS of the path length between the tips;  or 'mds_dist' to order the tips based on MDS of the distance between the tips (default: 'mds_dist')

+##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; 'mode' to order the tips by the most common order; 'mds' to order the tips based on MDS of the path length between the tips;  or 'mds_dist' to order the tips based on MDS of the distance between the tips (default: 'mode')

 ##' @param align.tips TRUE to align trees by their tips and FALSE to align trees by their root (default: TRUE)

 ##' @param jitter deviation to jitter tips

 ##' @param ... additional parameters passed to fortify, ggtree and geom_tree

@@ -16,41 +16,42 @@

 ##' @examples

 ##' require(ape)

 ##' require(dplyr)

+##' require(tidyr)

 ##' 

-##' # Plot mutliple trees with aligned tips

+##' # Plot multiple trees with aligned tips

 ##' trees <- list(read.tree(text="((a:1,b:1):1.5,c:2.5);"), read.tree(text="((a:1,c:1):1,b:2);"));

 ##' ggdensitree(trees) + geom_tiplab()

 ##' 

-##' # Plot multiple trees with aligmned tips with tip labls and separate tree colors

+##' # Plot multiple trees with aligned tips with tip labels and separate tree colors

 ##' trees.fort <- list(trees[[1]] %>% fortify %>% mutate(tree="a"), trees[[2]] %>% fortify %>% mutate(tree="b"));

 ##' ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black')

 ##' 

 ##' 

 ##' # Generate example data

 ##' set.seed(1)

-##' trees <- rmtree(5, 10)

-##' time.trees <- lapply(1:length(trees), function(i) {

-##'  	tree <- trees[[i]]

+##' random.trees <- rmtree(5, 10)

+##' time.trees <- lapply(seq_along(random.trees), function(i) {

+##'  	tree <- random.trees[[i]]

 ##'  	tree$tip.label <- paste0("t", 1:10)

 ##' 	dates <- estimate.dates(tree, 1:10, mu=1, nsteps=1)

 ##' 	tree$edge.length <- dates[tree$edge[, 2]] - dates[tree$edge[, 1]]

 ##' 	fortify(tree) %>% mutate(tree=factor(i, levels=as.character(1:10)))

 ##' })

 ##' 

-##' # Plot multiple trees with aligned tips from muliple time points

+##' # Plot multiple trees with aligned tips from multiple time points

 ##' ggdensitree(time.trees, aes(colour=tree), tip.order=paste0("t", 1:10)) + geom_tiplab(colour='black')

 ##' 

 ##' 

 ##' # Read example data

-##' trees <- read.tree(system.file("examples", "ggdensitree_example.tree", package="ggtree"))

+##' example.trees <- read.tree(system.file("examples", "ggdensitree_example.tree", package="ggtree"))

 ##' 

 ##' # Compute OTU

 ##' grp <- list(A = c("a.t1", "a.t2", "a.t3", "a.t4"), B = c("b.t1", "b.t2", "b.t3", "b.t4"), C = c("c.t1", "c.t2", "c.t3", "c.t4"))

-##' trees <- lapply(trees, groupOTU, grp)

+##' otu.trees <- lapply(example.trees, groupOTU, grp)

 ##' 

 ##' # Plot multiple trees colored by OTU

-##' ggdensitree(trees, aes(colour=group), alpha=1/6) + scale_colour_manual(values=c("black", "red", "green", "blue"))

-ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds_dist',

+##' ggdensitree(otu.trees, aes(colour=group), alpha=1/6, tip.order='mds') + scale_colour_manual(values=c("black", "red", "green", "blue"))

+ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mode',

 						align.tips=TRUE, jitter=0, ...) {

 	## factorize to simplify code

 	trees <- lapply(data, as.phylo)

@@ -59,7 +60,31 @@ ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='md

 	## determine tip order

 	if (length(tip.order) == 1) {

-		if (grepl('mds', tip.order)) {

+	  if (tip.order == 'mode') {

+	    first.label <- trees.f[[1]] %>%

+	      dplyr::filter(isTip) %>%

+	      dplyr::arrange(y) %>%

+	      dplyr::pull(label)

+

+		  res <- lapply(

+	      trees.f,

+	      . %>%

+	        dplyr::filter(isTip) %>%

+	        dplyr::arrange(y) %>%

+	        `$`("label") %>%

+	        match(first.label)

+	    ) %>%

+	      do.call(rbind, .) %>%

+	      as.data.frame() %>%

+	      dplyr::group_by_all() %>%

+	      dplyr::summarize(count = dplyr::n(), .groups = 'drop') %>%

+	      dplyr::filter(count == max(count)) %>%

+	      dplyr::select(-count) %>%

+	      dplyr::slice(1) %>%

+	      unlist()

+	    

+	    tip.order <- first.label[res]

+		} else if (grepl('mds', tip.order)) {

 			method <- tip.order

 			first.label <- trees.f[[1]] %>%

@@ -4,7 +4,7 @@

 ##' @param data a list of phylo objects or any object with an as.phylo and fortify method

 ##' @param mapping aesthetic mapping

 ##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')

-##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; or 'mds' to order the tips based on MDS of the distance between the tips (default: 'mds')

+##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; 'mds' to order the tips based on MDS of the path length between the tips;  or 'mds_dist' to order the tips based on MDS of the distance between the tips (default: 'mds')

 ##' @param align.tips TRUE to align trees by their tips and FALSE to align trees by their root (default: TRUE)

 ##' @param jitter deviation to jitter tips

 ##' @param ... additional parameters passed to fortify, ggtree and geom_tree

@@ -17,12 +17,16 @@

 ##' require(ape)

 ##' require(dplyr)

 ##' 

+##' # Plot mutliple trees with aligned tips

 ##' trees <- list(read.tree(text="((a:1,b:1):1.5,c:2.5);"), read.tree(text="((a:1,c:1):1,b:2);"));

 ##' ggdensitree(trees) + geom_tiplab()

 ##' 

+##' # Plot multiple trees with aligmned tips with tip labls and separate tree colors

 ##' trees.fort <- list(trees[[1]] %>% fortify %>% mutate(tree="a"), trees[[2]] %>% fortify %>% mutate(tree="b"));

 ##' ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black')

 ##' 

+##' 

+##' # Generate example data

 ##' set.seed(1)

 ##' trees <- rmtree(5, 10)

 ##' time.trees <- lapply(1:length(trees), function(i) {

@@ -31,9 +35,22 @@

 ##' 	dates <- estimate.dates(tree, 1:10, mu=1, nsteps=1)

 ##' 	tree$edge.length <- dates[tree$edge[, 2]] - dates[tree$edge[, 1]]

 ##' 	fortify(tree) %>% mutate(tree=factor(i, levels=as.character(1:10)))

-##' }) 

+##' })

+##' 

+##' # Plot multiple trees with aligned tips from muliple time points

 ##' ggdensitree(time.trees, aes(colour=tree), tip.order=paste0("t", 1:10)) + geom_tiplab(colour='black')

-ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds',

+##' 

+##' 

+##' # Read example data

+##' trees <- read.tree(system.file("ggtree", "inst/examples/ggdensitree_example.tree"))

+##' 

+##' # Compute OTU

+##' grp <- list(A = c("a.t1", "a.t2", "a.t3", "a.t4"), B = c("b.t1", "b.t2", "b.t3", "b.t4"), C = c("c.t1", "c.t2", "c.t3", "c.t4"))

+##' trees <- lapply(trees, groupOTU, grp)

+##' 

+##' # Plot multiple trees colored by OTU

+##' ggdensitree(trees, aes(colour=group), alpha=1/6) + scale_colour_manual(values=c("black", "red", "green", "blue"))

+ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds_dist',

 						align.tips=TRUE, jitter=0, ...) {

 	## factorize to simplify code

 	trees <- lapply(data, as.phylo)

@@ -42,7 +59,7 @@ ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='md

 	## determine tip order

 	if (length(tip.order) == 1) {

-		if (tip.order == 'mds') {

+		if (grepl('mds', tip.order)) {

 			first.label <- trees.f[[1]] %>%

 				dplyr::filter(.data$isTip) %>%

 				dplyr::pull(.data$label)

@@ -96,7 +113,7 @@ ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='md

 ## wrapper for cohpenetic to ensure that branch lengths exist

 cophenetic.phylo.check.length <- function(tree) {

-	if (is.null(tree$edge.length))

+	if (method != 'mds_dist' || is.null(tree$edge.length))

 		tree$edge.length <- rep(1, nrow(tree$edge))

 	stats::cophenetic(tree)

@@ -4,9 +4,8 @@

 ##' @param data a list of phylo objects or any object with an as.phylo and fortify method

 ##' @param mapping aesthetic mapping

 ##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')

-##' @param branch.length variable to be used to scale branch length. Setting to 'branch.length' will use the branch lengths of the tree objects. Default is 'none' to discard branch length and only plot cladogram (more reasonable for densitree).

 ##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; or 'mds' to order the tips based on MDS of the distance between the tips (default: 'mds')

-##' @param dangle TRUE to align trees by their tips and FALSE to align treesby their root (default: TRUE)

+##' @param align.tips TRUE to align trees by their tips and FALSE to align trees by their root (default: TRUE)

 ##' @param jitter deviation to jitter tips

 ##' @param ... additional parameters passed to fortify, ggtree and geom_tree

 ##' @return tree layer

@@ -17,78 +16,88 @@

 ##' @examples

 ##' require(ape)

 ##' require(dplyr)

-##' library(ape)

 ##' 

-##' trees <- list(read.tree(text="((a,b),c);"), read.tree(text="((a,c),b);"));

+##' trees <- list(read.tree(text="((a:1,b:1):1.5,c:2.5);"), read.tree(text="((a:1,c:1):1,b:2);"));

 ##' ggdensitree(trees) + geom_tiplab()

 ##' 

 ##' trees.fort <- list(trees[[1]] %>% fortify %>% mutate(tree="a"), trees[[2]] %>% fortify %>% mutate(tree="b"));

 ##' ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black')

-ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", branch.length = "none",

-                        tip.order='mds', dangle=TRUE, jitter=0, ...) {

-    ## factorize to simplify code

-    trees <- lapply(data, as.phylo)

-    trees.f <- lapply(data, fortify, layout=layout, branch.length = branch.length, ...)

-    n.tips <- sum(trees.f[[1]]$isTip)

-

-    ## determine tip order

-    if (length(tip.order) == 1) {

-        if (tip.order == 'mds') {

-            first.label <- trees.f[[1]] %>%

-                dplyr::filter(.data$isTip) %>%

-                dplyr::pull(.data$label)

-            

-            tip.order <- lapply(trees, function(x) {

-                match(x$tip.label, first.label)

-            })

-            

-            tip.2.tip <- lapply(trees, cophenetic.phylo.check.length)

-            tip.2.tip <- lapply(1:length(trees), function(i) {

-                tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]]

-            })

-            

-            all.tab <- do.call(rbind, tip.2.tip)

-            rownames(all.tab) <- NULL

-            

-            distances <- stats::dist(t(all.tab))

-            

-            res <- stats::cmdscale(distances, k=1)

+##' 

+##' set.seed(1)

+##' trees <- rmtree(5, 10)

+##' time.trees <- lapply(1:length(trees), function(i) {

+##'  	tree <- trees[[i]]

+##'  	tree$tip.label <- paste0("t", 1:10)

+##' 	dates <- estimate.dates(tree, 1:10, mu=1, nsteps=1)

+##' 	tree$edge.length <- dates[tree$edge[, 2]] - dates[tree$edge[, 1]]

+##' 	fortify(tree) %>% mutate(tree=factor(i, levels=as.character(1:10)))

+##' }) 

+##' ggdensitree(time.trees, aes(colour=tree), tip.order=paste0("t", 1:10)) + geom_tiplab(colour='black')

+ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds',

+						align.tips=TRUE, jitter=0, ...) {

+	## factorize to simplify code

+	trees <- lapply(data, as.phylo)

+	trees.f <- lapply(data, fortify, layout=layout, ...)

+	n.tips <- sum(trees.f[[1]]$isTip)

+	

+	## determine tip order

+	if (length(tip.order) == 1) {

+		if (tip.order == 'mds') {

+			first.label <- trees.f[[1]] %>%

+				dplyr::filter(.data$isTip) %>%

+				dplyr::pull(.data$label)

+			

+			tip.order <- lapply(trees, function(x) {

+				match(x$tip.label, first.label)

+			})

+			

+			tip.2.tip <- lapply(trees, cophenetic.phylo.check.length)

+			tip.2.tip <- lapply(1:length(trees), function(i) {

+				tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]]

+			})

+			

+			all.tab <- do.call(rbind, tip.2.tip)

+			rownames(all.tab) <- NULL

-            tip.order <- first.label[order(res[,1])]

-        } else if (as.numeric(tip.order) && tip.order <= length(trees)) {

-            labels <- trees.f[[tip.order]] %>%

-                dplyr::filter(.data$isTip) %>%

-                dplyr::pull(.data$label)

-            

-            tip.order <- labels[as.integer(trees.f[[tip.order]]$y)]

-        }

-    }

+			distances <- stats::dist(t(all.tab))

+			

+			res <- stats::cmdscale(distances, k=1)

+			

+			tip.order <- first.label[order(res[,1])]

+		} else if (as.numeric(tip.order) && tip.order <= length(trees)) {

+			labels <- trees.f[[tip.order]] %>%

+				dplyr::filter(.data$isTip) %>%

+				dplyr::pull(.data$label)

+			

+			tip.order <- labels[as.integer(trees.f[[tip.order]]$y)]

+		}

+	}

-    ## reorder tips (and shift x id dangling)

-    max.x <- vapply(trees.f, function(x) max(x$x, na.rm = TRUE), numeric(1))

-    farthest <- max(max.x)

-    trees.f <- lapply(1:length(trees), function(i) {

-        trees.f[[i]]$y <- getYcoord(trees[[i]], tip.order = tip.order)

-        if (dangle) {

-            trees.f[[i]]$x <- trees.f[[i]]$x + (farthest - max.x[i])

-        }

-        if (i > 1 && jitter > 0) {

-            trees.f[[i]]$y[1:n.tips] %<>% add(stats::rnorm(n.tips, mean=0, sd=jitter))

-        }

-        trees.f[[i]]

-    })

-    

-    ## plot all trees together

-    p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...)

-    for (x in trees.f[-1])

-        p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...)

-    p

+	## reorder tips (and shift x id align tips)

+	max.x <- vapply(trees.f, function(x) max(x$x, na.rm = TRUE), numeric(1))

+	farthest <- max(max.x)

+	trees.f <- lapply(1:length(trees), function(i) {

+		trees.f[[i]]$y <- getYcoord(trees[[i]], tip.order = tip.order)

+		if (align.tips) {

+			trees.f[[i]]$x <- trees.f[[i]]$x + (farthest - max.x[i])

+		}

+		if (i > 1 && jitter > 0) {

+			trees.f[[i]]$y[1:n.tips] %<>% add(stats::rnorm(n.tips, mean=0, sd=jitter))

+		}

+		trees.f[[i]]

+	})

+	

+	## plot all trees together

+	p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...)

+	for (x in trees.f[-1])

+		p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...)

+	p

 }

 ## wrapper for cohpenetic to ensure that branch lengths exist

 cophenetic.phylo.check.length <- function(tree) {

-    if (is.null(tree$edge.length))

-        tree$edge.length <- rep(1, nrow(tree$edge))

-    

-    stats::cophenetic(tree)

-}

+	if (is.null(tree$edge.length))

+		tree$edge.length <- rep(1, nrow(tree$edge))

+	

+	stats::cophenetic(tree)

+}

\ No newline at end of file

@@ -4,6 +4,7 @@

 ##' @param data a list of phylo objects or any object with an as.phylo and fortify method

 ##' @param mapping aesthetic mapping

 ##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')

+##' @param branch.length variable to be used to scale branch length. Setting to 'branch.length' will use the branch lengths of the tree objects. Default is 'none' to discard branch length and only plot cladogram (more reasonable for densitree).

 ##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; or 'mds' to order the tips based on MDS of the distance between the tips (default: 'mds')

 ##' @param dangle TRUE to align trees by their tips and FALSE to align treesby their root (default: TRUE)

 ##' @param jitter deviation to jitter tips

@@ -17,106 +18,77 @@

 ##' require(ape)

 ##' require(dplyr)

 ##' library(ape)

-##' library(tidyverse)

 ##' 

 ##' trees <- list(read.tree(text="((a,b),c);"), read.tree(text="((a,c),b);"));

 ##' ggdensitree(trees) + geom_tiplab()

 ##' 

 ##' trees.fort <- list(trees[[1]] %>% fortify %>% mutate(tree="a"), trees[[2]] %>% fortify %>% mutate(tree="b"));

 ##' ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black')

-ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds', dangle=TRUE, jitter=0, ...) {

-	# factorize to simplify code

-	trees <- lapply(data, as.phylo)

-	trees.f <- lapply(data, fortify, layout=layout, ...)

-	n.tips <- sum(trees.f[[1]]$isTip)

-	

-	# determine tip order

-	if (length(tip.order) == 1) {

-		if (tip.order == 'mds') {

-			first.label <- subset(trees.f[[1]], isTip)$label

-			

-			tip.order <- lapply(trees, . %$% tip.label %>% match(first.label))

-			

-			tip.2.tip <- lapply(trees, cophenetic.phylo.check.length)

-			tip.2.tip <- lapply(1:length(trees), function(i) tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]])

-			

-			all.tab <- do.call(rbind, tip.2.tip)

-			rownames(all.tab) <- NULL

-			

-			distances <- dist(t(all.tab))

-			

-			res <- cmdscale(distances, k=1)

-			

-			tip.order <- first.label[order(res[,1])]

-		} else if (as.numeric(tip.order) && tip.order <= length(trees)) {

-			labels <- subset(trees.f[[tip.order]], isTip)$label

+ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", branch.length = "none",

+                        tip.order='mds', dangle=TRUE, jitter=0, ...) {

+    ## factorize to simplify code

+    trees <- lapply(data, as.phylo)

+    trees.f <- lapply(data, fortify, layout=layout, branch.length = branch.length, ...)

+    n.tips <- sum(trees.f[[1]]$isTip)

+

+    ## determine tip order

+    if (length(tip.order) == 1) {

+        if (tip.order == 'mds') {

+            first.label <- trees.f[[1]] %>%

+                dplyr::filter(.data$isTip) %>%

+                dplyr::pull(.data$label)

+            

+            tip.order <- lapply(trees, function(x) {

+                match(x$tip.label, first.label)

+            })

+            

+            tip.2.tip <- lapply(trees, cophenetic.phylo.check.length)

+            tip.2.tip <- lapply(1:length(trees), function(i) {

+                tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]]

+            })

+            

+            all.tab <- do.call(rbind, tip.2.tip)

+            rownames(all.tab) <- NULL

+            

+            distances <- stats::dist(t(all.tab))

+            

+            res <- stats::cmdscale(distances, k=1)

-			tip.order <- labels[as.integer(trees.f[[tip.order]]$y)]

-		}

-	}

+            tip.order <- first.label[order(res[,1])]

+        } else if (as.numeric(tip.order) && tip.order <= length(trees)) {

+            labels <- trees.f[[tip.order]] %>%

+                dplyr::filter(.data$isTip) %>%

+                dplyr::pull(.data$label)

+            

+            tip.order <- labels[as.integer(trees.f[[tip.order]]$y)]

+        }

+    }

-	# reorder tips (and shift x id dangling)

+    ## reorder tips (and shift x id dangling)

+    max.x <- vapply(trees.f, function(x) max(x$x, na.rm = TRUE), numeric(1))

+    farthest <- max(max.x)

     trees.f <- lapply(1:length(trees), function(i) {

-    	trees.f[[i]]$y <- getYcoord_order(trees[[i]], tip.order)

-    	if (i > 1) {

-    		trees.f[[i]]$y[1:n.tips] %<>% add(rnorm(n.tips, mean=0, sd=jitter))

-    		if (dangle)

-   	    		trees.f[[i]]$x <- trees.f[[1]]$x[1] - trees.f[[i]]$x[match(trees.f[[1]]$label, trees.f[[i]]$label)][1] + trees.f[[i]]$x

-    	}

-    	trees.f[[i]]

+        trees.f[[i]]$y <- getYcoord(trees[[i]], tip.order = tip.order)

+        if (dangle) {

+            trees.f[[i]]$x <- trees.f[[i]]$x + (farthest - max.x[i])

+        }

+        if (i > 1 && jitter > 0) {

+            trees.f[[i]]$y[1:n.tips] %<>% add(stats::rnorm(n.tips, mean=0, sd=jitter))

+        }

+        trees.f[[i]]

     })

-	

-    # plot all trees together

-	p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...)

-	for (x in trees.f[-1])

-		p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...)

-	p

+    

+    ## plot all trees together

+    p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...)

+    for (x in trees.f[-1])

+        p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...)

+    p

 }

 ## wrapper for cohpenetic to ensure that branch lengths exist

 cophenetic.phylo.check.length <- function(tree) {

-	if (is.null(tree$edge.length))

-		tree$edge.length <- rep(1, nrow(tree$edge))

-	

-	cophenetic(tree)

-}

-

-## this is an adaptation of old code in ggtree

-## 

-##' @importFrom magrittr %>%

-##' @importFrom magrittr equals

-getYcoord_order <- function(tr, tip.order) {

-    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)

-    y[1:Ntip] <- match(tr$tip.label, tip.order)

-    y[-(1:Ntip)] <- NA

-

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

-    pvec[child] = parent

-

-    currentNode <- 1:Ntip

-    while(anyNA(y)) {

-        pNode <- unique(pvec[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)

-        })

-

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

-    }

-

-    return(y)

+    if (is.null(tree$edge.length))

+        tree$edge.length <- rep(1, nrow(tree$edge))

+    

+    stats::cophenetic(tree)

 }

new file mode 100644

@@ -0,0 +1,122 @@

+##' drawing phylogenetic trees from list of phylo objects

+##'

+##' @title ggdensitree

+##' @param data a list of phylo objects or any object with an as.phylo and fortify method

+##' @param mapping aesthetic mapping

+##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')

+##' @param tip.order the order of the tips by a character vector of taxa names; or an integer, N, to order the tips by the order of the tips in the Nth tree; or 'mds' to order the tips based on MDS of the distance between the tips (default: 'mds')

+##' @param dangle TRUE to align trees by their tips and FALSE to align treesby their root (default: TRUE)

+##' @param jitter deviation to jitter tips

+##' @param ... additional parameters passed to fortify, ggtree and geom_tree

+##' @return tree layer

+##' @importFrom magrittr %<>%

+##' @importFrom magrittr add

+##' @export

+##' @author Yu Guangchuang, Bradley R. Jones

+##' @examples

+##' require(ape)

+##' require(dplyr)

+##' library(ape)

+##' library(tidyverse)

+##' 

+##' trees <- list(read.tree(text="((a,b),c);"), read.tree(text="((a,c),b);"));

+##' ggdensitree(trees) + geom_tiplab()

+##' 

+##' trees.fort <- list(trees[[1]] %>% fortify %>% mutate(tree="a"), trees[[2]] %>% fortify %>% mutate(tree="b"));

+##' ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black')

+ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds', dangle=TRUE, jitter=0, ...) {

+	# factorize to simplify code

+	trees <- lapply(data, as.phylo)

+	trees.f <- lapply(data, fortify, layout=layout, ...)

+	n.tips <- sum(trees.f[[1]]$isTip)

+	

+	# determine tip order

+	if (length(tip.order) == 1) {

+		if (tip.order == 'mds') {

+			first.label <- subset(trees.f[[1]], isTip)$label

+			

+			tip.order <- lapply(trees, . %$% tip.label %>% match(first.label))

+			

+			tip.2.tip <- lapply(trees, cophenetic.phylo.check.length)

+			tip.2.tip <- lapply(1:length(trees), function(i) tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]])

+			

+			all.tab <- do.call(rbind, tip.2.tip)

+			rownames(all.tab) <- NULL

+			

+			distances <- dist(t(all.tab))

+			

+			res <- cmdscale(distances, k=1)

+			

+			tip.order <- first.label[order(res[,1])]

+		} else if (as.numeric(tip.order) && tip.order <= length(trees)) {

+			labels <- subset(trees.f[[tip.order]], isTip)$label

+			

+			tip.order <- labels[as.integer(trees.f[[tip.order]]$y)]

+		}

+	}

+	

+	# reorder tips (and shift x id dangling)

+    trees.f <- lapply(1:length(trees), function(i) {

+    	trees.f[[i]]$y <- getYcoord_order(trees[[i]], tip.order)

+    	if (i > 1) {

+    		trees.f[[i]]$y[1:n.tips] %<>% add(rnorm(n.tips, mean=0, sd=jitter))

+    		if (dangle)

+   	    		trees.f[[i]]$x <- trees.f[[1]]$x[1] - trees.f[[i]]$x[match(trees.f[[1]]$label, trees.f[[i]]$label)][1] + trees.f[[i]]$x

+    	}

+    	trees.f[[i]]

+    })

+	

+    # plot all trees together

+	p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...)

+	for (x in trees.f[-1])

+		p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...)

+	p

+}

+

+## wrapper for cohpenetic to ensure that branch lengths exist

+cophenetic.phylo.check.length <- function(tree) {

+	if (is.null(tree$edge.length))

+		tree$edge.length <- rep(1, nrow(tree$edge))

+	

+	cophenetic(tree)

+}

+

+## this is an adaptation of old code in ggtree

+## 

+##' @importFrom magrittr %>%

+##' @importFrom magrittr equals

+getYcoord_order <- function(tr, tip.order) {

+    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)

+    y[1:Ntip] <- match(tr$tip.label, tip.order)

+    y[-(1:Ntip)] <- NA

+

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

+    pvec[child] = parent

+

+    currentNode <- 1:Ntip

+    while(anyNA(y)) {

+        pNode <- unique(pvec[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)

+        })

+

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

+    }

+

+    return(y)

+}