Xavier Prudent Xavier Prudent - 1 year ago 74
R Question

Draw on a phylogeny edge with R

I want to draw a symbol (cross) anywhere along an edge of a phylogeny with R.

Let us take that tree:



and plot it with a symbol at the hg19 tip:

tree = read.tree("tree.mod")
points( 0.172365, 1, col="red", pch=20)

that is easy for the tips: the 'x' value is the full branch length, and 'y' is 1,2,3...

but I don't see how to do it for inner nodes like for instance loxAfr3-triMan1. I have the 'x' but cannot find the 'y'...

Any input welcome!

Answer Source

OK. I feel like I must be missing an easier way, but much of those code seems to be tucked away inside C functions called by the package's plotting function. Nevertheless, here are some functions in R that should be able to extract the x and y positions for particular nodes.

getphylo_x <- function(tree, node) {
    if(is.character(node)) {
        node <- which(c(tree$tip.label, tree$node.label)==node)
    pi <- tree$edge[tree$edge[,2]==node, 1]
    if (length(pi)) {
        ei<-which(tree$edge[,1]==pi & tree$edge[,2]==node)
        tree$edge.length[ei] + Recall(tree, pi)
    } else {
        if(!is.null(tree$root.edge)) {
        } else {

getphylo_y <- function(tree, node) {
    if(is.character(node)) {
        node <- which(c(tree$tip.label, tree$node.label)==node)
    ci <- tree$edge[tree$edge[,1]==node, 2]
    if (length(ci)==2) {
        mean(c(Recall(tree, ci[1]), Recall(tree, ci[2])))
    } else if (length(ci)==0) {
        Ntip <- length(tree$tip.label)
        which(tree$edge[tree$edge[, 2] <= Ntip, 2] == node)
    } else {
        stop(paste("error", length(ci)))

To use them, you just pass in your tree and your node name. For example

node <- "loxAfr3-triMan1"
x <- getphylo_x(tree, node)
y <- getphylo_y(tree, node)

points(x,y,pch=18, col="red", cex=2)

which produces

enter image description here

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download