0

我想用 R 在系统发育边缘的任何地方画一个符号（十字）。

让我们拿那棵树：

树.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'...

欢迎任何输入！

4

3 回答 3

2

好的。我觉得我一定错过了一种更简单的方法，但是这些代码中的大部分似乎都隐藏在包的绘图函数调用的 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)

产生

在此处输入图像描述

于 2014-09-02T17:27:01.560 回答
1

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

plot(tree)
locator()

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

于 2014-10-20T17:37:47.790 回答
1

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

鉴于 OP 在他们的评论之一中链接的树：

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)
于 2016-11-04T19:30:25.457 回答