用 R 在系统发育边缘绘制答案

【问题标题】：Draw on a phylogeny edge with R用 R 在系统发育边缘绘制
【发布时间】：2014-09-02 13:47:10
【问题描述】：

我想用 R 在系统发育边缘的任意位置画一个符号（十字）。

让我们拿走那棵树：

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;

并在 hg19 尖端用符号绘制它：

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

提示很简单：'x' 值是完整的分支长度，'y' 是 1,2,3...

但我不知道如何为内部节点执行此操作，例如 loxAfr3-triMan1。我有'x'，但找不到'y'...

欢迎任何意见！

【问题讨论】：

你是怎么设计的？ plot(tree) 似乎不起作用；它给出了错误“树中有单个（非拆分）节点；您可能需要使用 collapse.singles()”。而且你不能只打电话给points()而不先绘制一些东西。
复制粘贴有些问题，请使用该树：idisk-srv1.mpi-cbg.de/~prudent/tree.mod
那么你到底想要什么？您可以在 Paint 中创建一个模型或其他东西来明确您要在哪里绘制点吗？
请看一下：idisk-srv1.mpi-cbg.de/~prudent/exampeTree.jpeg
谢谢，这很有帮助。

标签： r phylogeny

【解决方案1】：

好的。我觉得我一定错过了一种更简单的方法，但是这些代码中的大部分似乎都隐藏在包的绘图函数调用的 C 函数中。不过，这里有一些 R 中的函数应该能够提取特定节点的 x 和 y 位置。

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)))
    }
}

要使用它们，您只需传入您的树和节点名称。例如

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


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

产生

【讨论】：

太棒了！这很完美。事实上，我也只是提出了一个解决方案，但更丑陋的 $ long 。非常感谢和问候来自法国的密歇根！
我迟到了，但在我看来 node.height() 和 node.depth.edgelength() 可以做你想做的事，而不需要使用自定义函数。
@jeanlain 我不知道这些功能。如果您写一个答案来说明如何为给定节点名称提取这些值，那将很有帮助。
这真是太棒了！感谢分享代码！任何想法如何对循环系统发育做同样的事情？

【解决方案2】：

先生。 Flick 的回答似乎完全涵盖了这一点，但我想我会补充一点，在打开树形图的情况下运行定位器功能将允许您获取用户指定点的 x,y 坐标。对于像 x,y 比例的放置这样的事情很有用。

plot(tree)
locator()

它将返回在按下 Esc 之前单击的所有点的 x 和 y 坐标。

【讨论】：

【解决方案3】：

这是一个不需要自定义函数的解决方案。

鉴于 OP 在他们的一个 cmets 中链接的树：

tree <- read.tree()   #file name isn't specified so ape will read what you paste in the console
((A:0.2,(B:0.1,C:0.1):0.1):0.1,((E:0.1,F:0.1):0.1,D:0.2):0.1);


plot(tree)
#Some random internal node:
node <- sample(tree$Nnode, 1) + length(tree$tip.label)  #nodes are not named in this tree
#if you want a (named) tip:
tip <- which(tree$tip.label == "A")
points(node.depth.edgelength(tree)[c(node, tip)], node.height(tree)[c(node,tip)], pch=18, col="red", cex=2)

【讨论】：