Browse code

Merge pull request #255 from brj1/master

small ggdensitree modifications

Guangchuang Yu authored on 12/09/2019 00:40:48 • GitHub committed on 12/09/2019 00:40:48
Showing2 changed files

... ...
@@ -4,9 +4,8 @@
4 4
 ##' @param data a list of phylo objects or any object with an as.phylo and fortify method
5 5
 ##' @param mapping aesthetic mapping
6 6
 ##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')
7
-##' @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).
8 7
 ##' @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')
9
-##' @param dangle TRUE to align trees by their tips and FALSE to align treesby their root (default: TRUE)
8
+##' @param align.tips TRUE to align trees by their tips and FALSE to align trees by their root (default: TRUE)
10 9
 ##' @param jitter deviation to jitter tips
11 10
 ##' @param ... additional parameters passed to fortify, ggtree and geom_tree
12 11
 ##' @return tree layer
... ...
@@ -17,78 +16,88 @@
17 16
 ##' @examples
18 17
 ##' require(ape)
19 18
 ##' require(dplyr)
20
-##' library(ape)
21 19
 ##' 
22
-##' trees <- list(read.tree(text="((a,b),c);"), read.tree(text="((a,c),b);"));
20
+##' trees <- list(read.tree(text="((a:1,b:1):1.5,c:2.5);"), read.tree(text="((a:1,c:1):1,b:2);"));
23 21
 ##' ggdensitree(trees) + geom_tiplab()
24 22
 ##' 
25 23
 ##' trees.fort <- list(trees[[1]] %>% fortify %>% mutate(tree="a"), trees[[2]] %>% fortify %>% mutate(tree="b"));
26 24
 ##' ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black')
27
-ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", branch.length = "none",
28
-                        tip.order='mds', dangle=TRUE, jitter=0, ...) {
29
-    ## factorize to simplify code
30
-    trees <- lapply(data, as.phylo)
31
-    trees.f <- lapply(data, fortify, layout=layout, branch.length = branch.length, ...)
32
-    n.tips <- sum(trees.f[[1]]$isTip)
33
-
34
-    ## determine tip order
35
-    if (length(tip.order) == 1) {
36
-        if (tip.order == 'mds') {
37
-            first.label <- trees.f[[1]] %>%
38
-                dplyr::filter(.data$isTip) %>%
39
-                dplyr::pull(.data$label)
40
-            
41
-            tip.order <- lapply(trees, function(x) {
42
-                match(x$tip.label, first.label)
43
-            })
44
-            
45
-            tip.2.tip <- lapply(trees, cophenetic.phylo.check.length)
46
-            tip.2.tip <- lapply(1:length(trees), function(i) {
47
-                tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]]
48
-            })
49
-            
50
-            all.tab <- do.call(rbind, tip.2.tip)
51
-            rownames(all.tab) <- NULL
52
-            
53
-            distances <- stats::dist(t(all.tab))
54
-            
55
-            res <- stats::cmdscale(distances, k=1)
25
+##' 
26
+##' set.seed(1)
27
+##' trees <- rmtree(5, 10)
28
+##' time.trees <- lapply(1:length(trees), function(i) {
29
+##'  	tree <- trees[[i]]
30
+##'  	tree$tip.label <- paste0("t", 1:10)
31
+##' 	dates <- estimate.dates(tree, 1:10, mu=1, nsteps=1)
32
+##' 	tree$edge.length <- dates[tree$edge[, 2]] - dates[tree$edge[, 1]]
33
+##' 	fortify(tree) %>% mutate(tree=factor(i, levels=as.character(1:10)))
34
+##' }) 
35
+##' ggdensitree(time.trees, aes(colour=tree), tip.order=paste0("t", 1:10)) + geom_tiplab(colour='black')
36
+ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds',
37
+						align.tips=TRUE, jitter=0, ...) {
38
+	## factorize to simplify code
39
+	trees <- lapply(data, as.phylo)
40
+	trees.f <- lapply(data, fortify, layout=layout, ...)
41
+	n.tips <- sum(trees.f[[1]]$isTip)
42
+	
43
+	## determine tip order
44
+	if (length(tip.order) == 1) {
45
+		if (tip.order == 'mds') {
46
+			first.label <- trees.f[[1]] %>%
47
+				dplyr::filter(.data$isTip) %>%
48
+				dplyr::pull(.data$label)
49
+			
50
+			tip.order <- lapply(trees, function(x) {
51
+				match(x$tip.label, first.label)
52
+			})
53
+			
54
+			tip.2.tip <- lapply(trees, cophenetic.phylo.check.length)
55
+			tip.2.tip <- lapply(1:length(trees), function(i) {
56
+				tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]]
57
+			})
58
+			
59
+			all.tab <- do.call(rbind, tip.2.tip)
60
+			rownames(all.tab) <- NULL
56 61
 			
57
-            tip.order <- first.label[order(res[,1])]
58
-        } else if (as.numeric(tip.order) && tip.order <= length(trees)) {
59
-            labels <- trees.f[[tip.order]] %>%
60
-                dplyr::filter(.data$isTip) %>%
61
-                dplyr::pull(.data$label)
62
-            
63
-            tip.order <- labels[as.integer(trees.f[[tip.order]]$y)]
64
-        }
65
-    }
62
+			distances <- stats::dist(t(all.tab))
63
+			
64
+			res <- stats::cmdscale(distances, k=1)
65
+			
66
+			tip.order <- first.label[order(res[,1])]
67
+		} else if (as.numeric(tip.order) && tip.order <= length(trees)) {
68
+			labels <- trees.f[[tip.order]] %>%
69
+				dplyr::filter(.data$isTip) %>%
70
+				dplyr::pull(.data$label)
71
+			
72
+			tip.order <- labels[as.integer(trees.f[[tip.order]]$y)]
73
+		}
74
+	}
66 75
 	
67
-    ## reorder tips (and shift x id dangling)
68
-    max.x <- vapply(trees.f, function(x) max(x$x, na.rm = TRUE), numeric(1))
69
-    farthest <- max(max.x)
70
-    trees.f <- lapply(1:length(trees), function(i) {
71
-        trees.f[[i]]$y <- getYcoord(trees[[i]], tip.order = tip.order)
72
-        if (dangle) {
73
-            trees.f[[i]]$x <- trees.f[[i]]$x + (farthest - max.x[i])
74
-        }
75
-        if (i > 1 && jitter > 0) {
76
-            trees.f[[i]]$y[1:n.tips] %<>% add(stats::rnorm(n.tips, mean=0, sd=jitter))
77
-        }
78
-        trees.f[[i]]
79
-    })
80
-    
81
-    ## plot all trees together
82
-    p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...)
83
-    for (x in trees.f[-1])
84
-        p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...)
85
-    p
76
+	## reorder tips (and shift x id align tips)
77
+	max.x <- vapply(trees.f, function(x) max(x$x, na.rm = TRUE), numeric(1))
78
+	farthest <- max(max.x)
79
+	trees.f <- lapply(1:length(trees), function(i) {
80
+		trees.f[[i]]$y <- getYcoord(trees[[i]], tip.order = tip.order)
81
+		if (align.tips) {
82
+			trees.f[[i]]$x <- trees.f[[i]]$x + (farthest - max.x[i])
83
+		}
84
+		if (i > 1 && jitter > 0) {
85
+			trees.f[[i]]$y[1:n.tips] %<>% add(stats::rnorm(n.tips, mean=0, sd=jitter))
86
+		}
87
+		trees.f[[i]]
88
+	})
89
+	
90
+	## plot all trees together
91
+	p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...)
92
+	for (x in trees.f[-1])
93
+		p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...)
94
+	p
86 95
 }
87 96
 
88 97
 ## wrapper for cohpenetic to ensure that branch lengths exist
89 98
 cophenetic.phylo.check.length <- function(tree) {
90
-    if (is.null(tree$edge.length))
91
-        tree$edge.length <- rep(1, nrow(tree$edge))
92
-    
93
-    stats::cophenetic(tree)
94
-}
99
+	if (is.null(tree$edge.length))
100
+		tree$edge.length <- rep(1, nrow(tree$edge))
101
+	
102
+	stats::cophenetic(tree)
103
+}
95 104
\ No newline at end of file
... ...
@@ -5,8 +5,7 @@
5 5
 \title{ggdensitree}
6 6
 \usage{
7 7
 ggdensitree(data = NULL, mapping = NULL, layout = "slanted",
8
-  branch.length = "none", tip.order = "mds", dangle = TRUE,
9
-  jitter = 0, ...)
8
+  tip.order = "mds", align.tips = TRUE, jitter = 0, ...)
10 9
 }
11 10
 \arguments{
12 11
 \item{data}{a list of phylo objects or any object with an as.phylo and fortify method}
... ...
@@ -15,11 +14,9 @@ ggdensitree(data = NULL, mapping = NULL, layout = "slanted",
15 14
 
16 15
 \item{layout}{one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')}
17 16
 
18
-\item{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).}
19
-
20 17
 \item{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')}
21 18
 
22
-\item{dangle}{TRUE to align trees by their tips and FALSE to align treesby their root (default: TRUE)}
19
+\item{align.tips}{TRUE to align trees by their tips and FALSE to align trees by their root (default: TRUE)}
23 20
 
24 21
 \item{jitter}{deviation to jitter tips}
25 22
 
... ...
@@ -34,13 +31,23 @@ drawing phylogenetic trees from list of phylo objects
34 31
 \examples{
35 32
 require(ape)
36 33
 require(dplyr)
37
-library(ape)
38 34
 
39
-trees <- list(read.tree(text="((a,b),c);"), read.tree(text="((a,c),b);"));
35
+trees <- list(read.tree(text="((a:1,b:1):1.5,c:2.5);"), read.tree(text="((a:1,c:1):1,b:2);"));
40 36
 ggdensitree(trees) + geom_tiplab()
41 37
 
42 38
 trees.fort <- list(trees[[1]] \%>\% fortify \%>\% mutate(tree="a"), trees[[2]] \%>\% fortify \%>\% mutate(tree="b"));
43 39
 ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black')
40
+
41
+set.seed(1)
42
+trees <- rmtree(5, 10)
43
+time.trees <- lapply(1:length(trees), function(i) {
44
+ 	tree <- trees[[i]]
45
+ 	tree$tip.label <- paste0("t", 1:10)
46
+	dates <- estimate.dates(tree, 1:10, mu=1, nsteps=1)
47
+	tree$edge.length <- dates[tree$edge[, 2]] - dates[tree$edge[, 1]]
48
+	fortify(tree) \%>\% mutate(tree=factor(i, levels=as.character(1:10)))
49
+}) 
50
+ggdensitree(time.trees, aes(colour=tree), tip.order=paste0("t", 1:10)) + geom_tiplab(colour='black')
44 51
 }
45 52
 \author{
46 53
 Yu Guangchuang, Bradley R. Jones