|
| 1 | +#' @importFrom magrittr %>% |
| 2 | +#' @export |
| 3 | +magrittr::'%>%' |
| 4 | + |
| 5 | +#' @method fortify evonet |
| 6 | +#' @importFrom ggplot2 fortify |
| 7 | +#' @export |
| 8 | +fortify.evonet <- function(model, data, layout="rectangular", ladderize=FALSE, |
| 9 | + right=FALSE, mrsd=NULL, as.Date =FALSE, ...){ |
| 10 | + class(model) <- "phylo" |
| 11 | +# ggtree:::fortify.phylo |
| 12 | + df <- fortify(model, ladderize=ladderize) |
| 13 | + |
| 14 | + hybridEdge <- logical(nrow(df)) |
| 15 | + hybridEdge[grep("#", df$label)] <- TRUE |
| 16 | + df <- cbind(df, hybridEdge=hybridEdge) |
| 17 | + |
| 18 | + reticulation <- model$reticulation |
| 19 | + df.ret <- df[reticulation[,1], , drop=FALSE] |
| 20 | +# df.ret <- df[reticulation[,2], , drop=FALSE] |
| 21 | + df.ret[,c("node", "parent")] <- reticulation |
| 22 | + df.ret[, "hybridEdge"] <- TRUE |
| 23 | + df <- rbind(df, df.ret) |
| 24 | + df |
| 25 | +} |
| 26 | + |
| 27 | + |
| 28 | +#' drawing phylogenetic tree from phylo object |
| 29 | +#' |
| 30 | +#' |
| 31 | +#' @title ggevonet |
| 32 | +#' @param tr phylo object |
| 33 | +#' @param mapping aes mapping |
| 34 | +#' @param layout one of 'rectangular', 'slanted' |
| 35 | +#' @param open.angle open angle, only for 'fan' layout |
| 36 | +#' @param mrsd most recent sampling date |
| 37 | +#' @param as.Date logical whether using Date class in time tree |
| 38 | +#' @param yscale y scale |
| 39 | +#' @param yscale_mapping yscale mapping for category variable |
| 40 | +#' @param ladderize logical |
| 41 | +#' @param right logical |
| 42 | +#' @param branch.length variable for scaling branch, if 'none' draw cladogram |
| 43 | +#' @param ndigits number of digits to round numerical annotation variable |
| 44 | +#' @param min_crossing logical, rotate clades to minimize crossings |
| 45 | +#' @param ... additional parameter |
| 46 | +#' @return tree |
| 47 | +#' @seealso \code{\link[ape]{evonet}}, \code{\link[ggtree]{ggtree}} |
| 48 | +#' @importFrom ggplot2 ggplot |
| 49 | +#' @importFrom ggplot2 xlab |
| 50 | +#' @importFrom ggplot2 ylab |
| 51 | +#' @importFrom ggplot2 annotate |
| 52 | +#' @importFrom ggplot2 scale_x_reverse |
| 53 | +#' @importFrom ggplot2 ylim |
| 54 | +#' @importFrom ggplot2 coord_flip |
| 55 | +#' @importFrom ggplot2 coord_polar |
| 56 | +#' @importFrom ggplot2 aes |
| 57 | +#' @importFrom ggplot2 aes_ |
| 58 | +#' @importFrom ggtree geom_tree2 |
| 59 | +#' @importFrom ggtree theme_tree |
| 60 | +#' @importFrom phangorn coords |
| 61 | +#' @importFrom ape read.evonet |
| 62 | +#' @author Klaus Schliep |
| 63 | +#' @examples |
| 64 | +#' (enet <- ape::read.evonet(text="((a:2,(b:1)#H1:1):1,(#H1,c:1):2);")) |
| 65 | +#' ggevonet(enet) + geom_tiplab() |
| 66 | +#' @export |
| 67 | +ggevonet <- function (tr, mapping=NULL, layout="slanted", open.angle=0, |
| 68 | + mrsd=NULL, as.Date=FALSE, yscale="none", yscale_mapping=NULL, |
| 69 | + ladderize=FALSE, right=FALSE, branch.length="branch.length", |
| 70 | + ndigits=NULL, min_crossing=TRUE, ...) |
| 71 | +{ |
| 72 | + layout <- match.arg(layout, c("rectangular", "slanted")) |
| 73 | +# , "fan", "circular", "radial", "unrooted", "equal_angle", "daylight" |
| 74 | + |
| 75 | + if(is.null(tr$edge.length)){ |
| 76 | + nh <- node.depth.evonet(tr) |
| 77 | + tr$edge.length <- nh[tr$edge[,1]] - nh[tr$edge[,2]] |
| 78 | + } |
| 79 | + if(min_crossing){ |
| 80 | + tr <- minimize_overlap(tr) |
| 81 | + } |
| 82 | + if (yscale != "none") { |
| 83 | + layout <- "slanted" |
| 84 | + } |
| 85 | + if (is.null(mapping)) { |
| 86 | + mapping <- aes_(~x, ~y) |
| 87 | + } |
| 88 | + else { |
| 89 | + mapping <- modifyList(aes_(~x, ~y), mapping) |
| 90 | + } |
| 91 | + mapping <- modifyList(aes(linetype = hybridEdge), mapping) |
| 92 | + p <- ggplot(tr, mapping=mapping, layout=layout, mrsd=mrsd, as.Date=as.Date, |
| 93 | + yscale=yscale, yscale_mapping=yscale_mapping, |
| 94 | + ladderize=ladderize, right=right, branch.length=branch.length, |
| 95 | + ndigits=ndigits, ...) |
| 96 | + |
| 97 | + p <- p + geom_tree2(layout=layout, ...) |
| 98 | + p <- p + theme_tree(legend.position="none") |
| 99 | + class(p) <- c("ggtree", class(p)) |
| 100 | + return(p) |
| 101 | +} |
| 102 | + |
| 103 | + |
| 104 | +#' @title minimize_overlap |
| 105 | +#' reduces reticulation lines crossing over in plots |
| 106 | +#' @param x Tree of class 'evonet' |
| 107 | +#' @return A Tree with rotated nodes of class 'evonet' |
| 108 | +#' @importFrom ape node.height rotate Ntip |
| 109 | +#' @importFrom phangorn Ancestors |
| 110 | +#' @author L. Francisco Henao Diaz |
| 111 | +#' @examples |
| 112 | +#' fishnet <- ape::read.evonet(text="(Xalvarezi,Xmayae,((Xsignum,((Xmonticolus, |
| 113 | +#' (Xclemenciae_F2,#H25)),(((((((((Xgordoni,Xmeyeri),Xcouchianus),Xvariatus), |
| 114 | +#' Xevelynae),(Xxiphidium,#H24)),Xmilleri),Xandersi),Xmaculatus),(((Xmontezumae, |
| 115 | +#' (Xcortezi,(Xbirchmanni_GARC,Xmalinche_CHIC2))),((Xnigrensis,Xmultilineatus), |
| 116 | +#' (Xpygmaeus,Xcontinens))))#H24))),(Xhellerii)#H25));") |
| 117 | +#' fishnet$edge.length <- NULL |
| 118 | +#' new_tre <- minimize_overlap(fishnet) |
| 119 | +#' |
| 120 | +#' par(mfrow=c(1,2)) |
| 121 | +#' ggevonet(fishnet) |
| 122 | +#' ggevonet(new_tre) |
| 123 | +#' |
| 124 | +#' net2 <- ape::read.evonet(text="(15,(1,((14,(#H1,(((12,13),(11,#H3)),(7, |
| 125 | +#' ((10)#H3,(8,9)))))),((((2,3))#H2,(6,(5,(#H2,4)))))#H1)));") |
| 126 | +#' # Cui et al. 2013 Evol. |
| 127 | +#' new_net2 <- minimize_overlap(net2) |
| 128 | +#' ggevonet(net2) |
| 129 | +#' ggevonet(new_net2) |
| 130 | +#' |
| 131 | +#' @export |
| 132 | +minimize_overlap <- function(x){ |
| 133 | + if(class(x)[1]!='evonet') stop("x should be an 'evonet' class") |
| 134 | + n_iter <- round(x$Nnode*3/4) |
| 135 | + # r_hist <- c() |
| 136 | + for(j in seq_len(n_iter)){ |
| 137 | + h <- node.height(x) |
| 138 | + best_r <- sum(abs(h[x$reticulation[,1]]- h[x$reticulation[,2]])) |
| 139 | + best_c <- -1 |
| 140 | + # r_hist[j] <- best_r |
| 141 | + nodes2rot <- intersect(sort(unique(unlist(Ancestors(x, |
| 142 | + c(x$reticulation))))), which(tabulate(x$edge[,1]) > 1) ) |
| 143 | + for(i in seq_along(nodes2rot)){ |
| 144 | + nh <- node.height(ape::rotate(x, nodes2rot[i])) |
| 145 | + best_nr <- sum(abs(nh[x$reticulation[,1]] - nh[x$reticulation[,2]])) |
| 146 | + if(best_nr < best_r){ |
| 147 | + best_c <- nodes2rot[i] |
| 148 | + best_r <- best_nr |
| 149 | + } |
| 150 | + } |
| 151 | + if(best_c > 0) x <- ape::rotate(x, best_c) |
| 152 | + else break() |
| 153 | + } |
| 154 | + x |
| 155 | +} |
| 156 | + |
| 157 | +#' These functions return the depths or heights of nodes and tips. |
| 158 | +#' |
| 159 | +#' @title Depth of Nodes |
| 160 | +#' @param x an object of class "evonet" |
| 161 | +#' @param \dots Further arguments passed to or from other methods. |
| 162 | +#' @return a vector with the depth of the nodes |
| 163 | +#' @seealso \code{\link[ape]{node.depth}} |
| 164 | +#' @examples |
| 165 | +#' z <- ape::read.evonet(text = "((1,((2,(3,(4)Y#H1)g)e, |
| 166 | +#' (((Y#H1, 5)h,6)f)X#H2)c)a,((X#H2,7)d,8)b)r;") |
| 167 | +#' nd <- node.depth.evonet(z) |
| 168 | +#' z$edge.length <- nd[z$edge[,1]] - nd[z$edge[,2]] |
| 169 | +#' ggevonet(z) |
| 170 | +#' |
| 171 | +#' @importFrom phangorn getRoot Ancestors Descendants |
| 172 | +#' @export |
| 173 | +node.depth.evonet <- function(x, ...){ |
| 174 | + x <- ape::reorder.phylo(x) |
| 175 | + root <- getRoot(x) |
| 176 | + max_nodes <- max(x$edge) |
| 177 | + nTip <- Ntip(x) |
| 178 | + desc <- Descendants(x, seq_len(max_nodes), "children") |
| 179 | + anc <- Ancestors(x) |
| 180 | + pa <- vector("list", max_nodes) |
| 181 | + ind <- which(x$edge[,2] > Ntip(x)) |
| 182 | + pa[x$edge[ind,2]] <- x$edge[ind,1] |
| 183 | + for(i in seq_len(nrow(x$reticulation))){ |
| 184 | + pa[[x$reticulation[i,2]]] <- sort( c(pa[[x$reticulation[i,2]]], |
| 185 | + x$reticulation[i,1] ) ) |
| 186 | + # pa[[x$reticulation[i,1] ]] <- numeric(0) |
| 187 | + } |
| 188 | + ind <- which(lengths(pa) > 0) |
| 189 | + depth <- numeric(max_nodes) |
| 190 | + depth[root] <- 1 |
| 191 | + done <- logical(max_nodes) |
| 192 | + done[root] <- TRUE |
| 193 | + candidates <- desc[[root]] |
| 194 | + candidates <- candidates[candidates>nTip] |
| 195 | + d <- 1 |
| 196 | + while(length(candidates)>0){ |
| 197 | + active <- vapply(candidates, function(x) all(done[pa[[x]]]), FALSE) |
| 198 | + tmp <- which(active)[1] #sample(active,1) |
| 199 | + candidates <- c(candidates, desc[[candidates[tmp] ]]) |
| 200 | + candidates <- candidates[candidates>nTip] |
| 201 | + d <- d+1 |
| 202 | + done[candidates[tmp]] <- TRUE |
| 203 | + depth[candidates[tmp]] <- d |
| 204 | + candidates <- candidates[-tmp] |
| 205 | + } |
| 206 | + depth <- d+2 - depth |
| 207 | + depth[seq_len(nTip)] <- 1 |
| 208 | + depth |
| 209 | +} |
0 commit comments