Added ggdensitree method
... | ... |
@@ -88,6 +88,7 @@ export(get.path) |
88 | 88 |
export(get_clade_position) |
89 | 89 |
export(get_heatmap_column_position) |
90 | 90 |
export(get_taxa_name) |
91 |
+export(ggdensitree) |
|
91 | 92 |
export(ggplot) |
92 | 93 |
export(ggsave) |
93 | 94 |
export(ggtree) |
... | ... |
@@ -223,6 +224,7 @@ importFrom(grid,unit) |
223 | 224 |
importFrom(grid,viewport) |
224 | 225 |
importFrom(magrittr,"%<>%") |
225 | 226 |
importFrom(magrittr,"%>%") |
227 |
+importFrom(magrittr,add) |
|
226 | 228 |
importFrom(magrittr,equals) |
227 | 229 |
importFrom(methods,is) |
228 | 230 |
importFrom(methods,missingArg) |
229 | 231 |
new file mode 100644 |
... | ... |
@@ -0,0 +1,122 @@ |
1 |
+##' drawing phylogenetic trees from list of phylo objects |
|
2 |
+##' |
|
3 |
+##' @title ggdensitree |
|
4 |
+##' @param data a list of phylo objects or any object with an as.phylo and fortify method |
|
5 |
+##' @param mapping aesthetic mapping |
|
6 |
+##' @param layout one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted') |
|
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') |
|
8 |
+##' @param dangle TRUE to align trees by their tips and FALSE to align treesby their root (default: TRUE) |
|
9 |
+##' @param jitter deviation to jitter tips |
|
10 |
+##' @param ... additional parameters passed to fortify, ggtree and geom_tree |
|
11 |
+##' @return tree layer |
|
12 |
+##' @importFrom magrittr %<>% |
|
13 |
+##' @importFrom magrittr add |
|
14 |
+##' @export |
|
15 |
+##' @author Yu Guangchuang, Bradley R. Jones |
|
16 |
+##' @examples |
|
17 |
+##' require(ape) |
|
18 |
+##' require(dplyr) |
|
19 |
+##' library(ape) |
|
20 |
+##' library(tidyverse) |
|
21 |
+##' |
|
22 |
+##' trees <- list(read.tree(text="((a,b),c);"), read.tree(text="((a,c),b);")); |
|
23 |
+##' ggdensitree(trees) + geom_tiplab() |
|
24 |
+##' |
|
25 |
+##' trees.fort <- list(trees[[1]] %>% fortify %>% mutate(tree="a"), trees[[2]] %>% fortify %>% mutate(tree="b")); |
|
26 |
+##' ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black') |
|
27 |
+ggdensitree <- function(data=NULL, mapping=NULL, layout="slanted", tip.order='mds', dangle=TRUE, jitter=0, ...) { |
|
28 |
+ # factorize to simplify code |
|
29 |
+ trees <- lapply(data, as.phylo) |
|
30 |
+ trees.f <- lapply(data, fortify, layout=layout, ...) |
|
31 |
+ n.tips <- sum(trees.f[[1]]$isTip) |
|
32 |
+ |
|
33 |
+ # determine tip order |
|
34 |
+ if (length(tip.order) == 1) { |
|
35 |
+ if (tip.order == 'mds') { |
|
36 |
+ first.label <- subset(trees.f[[1]], isTip)$label |
|
37 |
+ |
|
38 |
+ tip.order <- lapply(trees, . %$% tip.label %>% match(first.label)) |
|
39 |
+ |
|
40 |
+ tip.2.tip <- lapply(trees, cophenetic.phylo.check.length) |
|
41 |
+ tip.2.tip <- lapply(1:length(trees), function(i) tip.2.tip[[i]][tip.order[[i]], tip.order[[i]]]) |
|
42 |
+ |
|
43 |
+ all.tab <- do.call(rbind, tip.2.tip) |
|
44 |
+ rownames(all.tab) <- NULL |
|
45 |
+ |
|
46 |
+ distances <- dist(t(all.tab)) |
|
47 |
+ |
|
48 |
+ res <- cmdscale(distances, k=1) |
|
49 |
+ |
|
50 |
+ tip.order <- first.label[order(res[,1])] |
|
51 |
+ } else if (as.numeric(tip.order) && tip.order <= length(trees)) { |
|
52 |
+ labels <- subset(trees.f[[tip.order]], isTip)$label |
|
53 |
+ |
|
54 |
+ tip.order <- labels[as.integer(trees.f[[tip.order]]$y)] |
|
55 |
+ } |
|
56 |
+ } |
|
57 |
+ |
|
58 |
+ # reorder tips (and shift x id dangling) |
|
59 |
+ trees.f <- lapply(1:length(trees), function(i) { |
|
60 |
+ trees.f[[i]]$y <- getYcoord_order(trees[[i]], tip.order) |
|
61 |
+ if (i > 1) { |
|
62 |
+ trees.f[[i]]$y[1:n.tips] %<>% add(rnorm(n.tips, mean=0, sd=jitter)) |
|
63 |
+ if (dangle) |
|
64 |
+ 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 |
|
65 |
+ } |
|
66 |
+ trees.f[[i]] |
|
67 |
+ }) |
|
68 |
+ |
|
69 |
+ # plot all trees together |
|
70 |
+ p <- ggtree(tr=trees.f[[1]], mapping=mapping, layout=layout, ...) |
|
71 |
+ for (x in trees.f[-1]) |
|
72 |
+ p <- p + geom_tree(mapping=mapping, data=x, layout=layout, ...) |
|
73 |
+ p |
|
74 |
+} |
|
75 |
+ |
|
76 |
+## wrapper for cohpenetic to ensure that branch lengths exist |
|
77 |
+cophenetic.phylo.check.length <- function(tree) { |
|
78 |
+ if (is.null(tree$edge.length)) |
|
79 |
+ tree$edge.length <- rep(1, nrow(tree$edge)) |
|
80 |
+ |
|
81 |
+ cophenetic(tree) |
|
82 |
+} |
|
83 |
+ |
|
84 |
+## this is an adaptation of old code in ggtree |
|
85 |
+## |
|
86 |
+##' @importFrom magrittr %>% |
|
87 |
+##' @importFrom magrittr equals |
|
88 |
+getYcoord_order <- function(tr, tip.order) { |
|
89 |
+ Ntip <- length(tr[["tip.label"]]) |
|
90 |
+ N <- getNodeNum(tr) |
|
91 |
+ |
|
92 |
+ edge <- tr[["edge"]] |
|
93 |
+ parent <- edge[,1] |
|
94 |
+ child <- edge[,2] |
|
95 |
+ |
|
96 |
+ cl <- split(child, parent) |
|
97 |
+ child_list <- list() |
|
98 |
+ child_list[as.numeric(names(cl))] <- cl |
|
99 |
+ |
|
100 |
+ y <- numeric(N) |
|
101 |
+ y[1:Ntip] <- match(tr$tip.label, tip.order) |
|
102 |
+ y[-(1:Ntip)] <- NA |
|
103 |
+ |
|
104 |
+ pvec <- integer(max(tr$edge)) |
|
105 |
+ pvec[child] = parent |
|
106 |
+ |
|
107 |
+ currentNode <- 1:Ntip |
|
108 |
+ while(anyNA(y)) { |
|
109 |
+ pNode <- unique(pvec[currentNode]) |
|
110 |
+ |
|
111 |
+ idx <- sapply(pNode, function(i) all(child_list[[i]] %in% currentNode)) |
|
112 |
+ newNode <- pNode[idx] |
|
113 |
+ |
|
114 |
+ y[newNode] <- sapply(newNode, function(i) { |
|
115 |
+ mean(y[child_list[[i]]], na.rm=TRUE) |
|
116 |
+ }) |
|
117 |
+ |
|
118 |
+ currentNode <- c(currentNode[!currentNode %in% unlist(child_list[newNode])], newNode) |
|
119 |
+ } |
|
120 |
+ |
|
121 |
+ return(y) |
|
122 |
+} |
0 | 123 |
new file mode 100644 |
... | ... |
@@ -0,0 +1,45 @@ |
1 |
+% Generated by roxygen2: do not edit by hand |
|
2 |
+% Please edit documentation in R/ggdensitree.R |
|
3 |
+\name{ggdensitree} |
|
4 |
+\alias{ggdensitree} |
|
5 |
+\title{ggdensitree} |
|
6 |
+\usage{ |
|
7 |
+ggdensitree(data = NULL, mapping = NULL, layout = "slanted", |
|
8 |
+ tip.order = "mds", dangle = TRUE, jitter = 0, ...) |
|
9 |
+} |
|
10 |
+\arguments{ |
|
11 |
+\item{data}{a list of phylo objects or any object with an as.phylo and fortify method} |
|
12 |
+ |
|
13 |
+\item{mapping}{aesthetic mapping} |
|
14 |
+ |
|
15 |
+\item{layout}{one of 'slanted', 'rectangluar', 'fan', 'circular' or 'radial' (default: 'slanted')} |
|
16 |
+ |
|
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')} |
|
18 |
+ |
|
19 |
+\item{dangle}{TRUE to align trees by their tips and FALSE to align treesby their root (default: TRUE)} |
|
20 |
+ |
|
21 |
+\item{jitter}{deviation to jitter tips} |
|
22 |
+ |
|
23 |
+\item{...}{additional parameters passed to fortify, ggtree and geom_tree} |
|
24 |
+} |
|
25 |
+\value{ |
|
26 |
+tree layer |
|
27 |
+} |
|
28 |
+\description{ |
|
29 |
+drawing phylogenetic trees from list of phylo objects |
|
30 |
+} |
|
31 |
+\examples{ |
|
32 |
+require(ape) |
|
33 |
+require(dplyr) |
|
34 |
+library(ape) |
|
35 |
+library(tidyverse) |
|
36 |
+ |
|
37 |
+trees <- list(read.tree(text="((a,b),c);"), read.tree(text="((a,c),b);")); |
|
38 |
+ggdensitree(trees) + geom_tiplab() |
|
39 |
+ |
|
40 |
+trees.fort <- list(trees[[1]] \%>\% fortify \%>\% mutate(tree="a"), trees[[2]] \%>\% fortify \%>\% mutate(tree="b")); |
|
41 |
+ggdensitree(trees.fort, aes(colour=tree)) + geom_tiplab(colour='black') |
|
42 |
+} |
|
43 |
+\author{ |
|
44 |
+Yu Guangchuang, Bradley R. Jones |
|
45 |
+} |