Draw on a phylogeny edge with R

Question

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

Let us take that tree:

tree.mod

(((((hg19:0.0295977,macFas5:0.0351997)hg19-macFas5:0.0796862,otoGar3:0.153768)hg19-otoGar3:0.0189364,(((speTri2:0.136745,(mm10:0.0836004,rn5:0.0894755)mm10-rn5:0.221692)speTri2-mm10:0.009583,cavPor3:0.230585)speTri2-cavPor3:0.0279417,oryCun2:0.212423)speTri2-oryCun2:0.0141334)hg19-speTri2:0.0224427,((susScr3:0.127502,(camFer1:0.111428((HLbalAcu1:0.0196864,(HLphyCat1:0.0216245,(HLlipVex1:0.0178214,(turTru2:0.00660779,orcOrc1:0.00557822)turTru2-orcOrc1:0.0113362)HLlipVex1-turTru2:0.00784223)HLphyCat1-HLlipVex1:0.00385294)HLbalAcu1-HLphyCat1:0.0417511(bosTau7:0.0320796,oviAri3:0.0370855)bosTau7-oviAri3:0.0963575)HLbalAcu1-bosTau7:0.0212413)camFer1-HLbalAcu1:0.00459573)susScr3-camFer1:0.0425276,((equCab2:0.107972,(felCat5:0.0893591,(canFam3:0.0908515,((ailMel1:0.0208923,HLursMar1:0.0194161)ailMel1-HLursMar1:0.041512,(odoRosDiv1:0.0259201,lepWed1:0.0238902)odoRosDiv1-lepWed1:0.0288135)ailMel1-odoRosDiv1:0.0213261)canFam3-ailMel1:0.0207311)felCat5-canFam3:0.0501969)equCab2-felCat5:0.00485202,myoLuc2:0.169981)equCab2-myoLuc2:0.00361502)susScr3-equCab2:0.0302374)hg19-susScr3:0.0217021,(((loxAfr3:0.0821245,triMan1:0.0662932)loxAfr3-triMan1:0.0296365,echTel2:0.245195)loxAfr3-echTel2:0.0492624,dasNov3:0.159632)loxAfr3-dasNov3:0.00825218)hg19-loxAfr3;

and plot it with a symbol at the hg19 tip:

R
library(ape)
tree = read.tree("tree.mod")
plot(tree)
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 1

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)) {
            tree$root.edge
        } else {
            0
        }
    }
}

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)


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

which produces